Novel domain expansion methods to improve the computational efficiency of the Chemical Master Equation solution for large biological networks

Background Numerical solutions of the chemical master equation (CME) are important for understanding the stochasticity of biochemical systems. However, solving CMEs is a formidable task. This task is complicated due to the nonlinear nature of the reactions and the size of the networks which result in different realizations. Most importantly, the exponential growth of the size of the state-space, with respect to the number of different species in the system makes this a challenging assignment. When the biochemical system has a large number of variables, the CME solution becomes intractable. We introduce the intelligent state projection (ISP) method to use in the stochastic analysis of these systems. For any biochemical reaction network, it is important to capture more than one moment: this allows one to describe the system’s dynamic behaviour. ISP is based on a state-space search and the data structure standards of artificial intelligence (AI). It can be used to explore and update the states of a biochemical system. To support the expansion in ISP, we also develop a Bayesian likelihood node projection (BLNP) function to predict the likelihood of the states. Results To demonstrate the acceptability and effectiveness of our method, we apply the ISP method to several biological models discussed in prior literature. The results of our computational experiments reveal that the ISP method is effective both in terms of the speed and accuracy of the expansion, and the accuracy of the solution. This method also provides a better understanding of the state-space of the system in terms of blueprint patterns. Conclusions The ISP is the de-novo method which addresses both accuracy and performance problems for CME solutions. It systematically expands the projection space based on predefined inputs. This ensures accuracy in the approximation and an exact analytical solution for the time of interest. The ISP was more effective both in predicting the behavior of the state-space of the system and in performance management, which is a vital step towards modeling large biochemical systems.


Background
In systems biology, it is crucial to understand the dynamics of large and complicated biochemical reaction networks. Recent advances in computing and mathematical techniques mean it is easier for biologists to deal with enormous amounts of experimental data, right down to the level of a single molecule of a species. Such information reveals the presence of a high level of stochasticity in the networks of biochemical reactions. In biochemical reaction networks, stochastic models have made significant contributions to the fields of systems biology [1,2], neuroscience [3], and drug modeling [4]. In a complex system, biochemical reactions are often modeled as reaction rate equations (RREs) using ordinary differential equations (ODEs). Examples of this kind of work include the biochemical networks of Alzheimer's disease (AD) [5]; the pathways in the fungal pathogen Candida albicans [6]; and the COVID-19 coronavirus pathogen network [7]. In each of these examples, the behavior of different pathways is still largely unknown. All these models only contain species with small copy numbers and widely different reaction rates; the probabilistic descriptions of time evolution of molecular concentrations (or numbers) are more suited for understanding the dynamics of such systems. One probabilistic approach for modeling a biochemical reaction network is to deduce a set of integro-differential equations known as chemical master equations (CMEs) [8,9]. CMEs describe the evolution of the probability distribution over the entire state-space of a biochemical system that jumps from one set of states to another set of states in continuous time: they are a continuous time version of Markov chains (CTMCs) [8,10] with discrete states. By defining the Markov chain [10,11], we can consider the joint and marginal probability densities of the species in a system that changes over time [12].
In such cases, the development of RREs with molecular numbers becomes very important. The biochemical reaction network can be defined in terms of the discrete state X ≡ ðx 1 ; …; xÑÞ T vector of non-negative integers xÑ for the given conditions, whereÑ ≥1. {X(t) : t ∈ K; φ} defines a stochastic process, where K is the indexing scheme and φ is the sample space. Following the derivation in [9], for every reaction, there exists a reaction channel, R M , which determines the unique reaction in the system with a propensity function k M . The specific combinations of the reactant species in R M will react during an infinitesimal [t, t + dt) time interval. The average probability a μ (X(t))dt of a particular R M fires within [t, t + dt) is the multiplication of the numbers of reactant species, denoted by square brackets, by k M . For example, In the case where the reactants are of the same type, for example A þ A → k 2 T , then a 2 :k 2 ð ½A½ðA − 1Þ 2 Þ. The set consisting of all the reaction channels, R M , is the union of sets of fast reactions and slow reactions [12]. They are categorized into sets of R M(fs) and R M(sr) reactions, respectively, based on their propensity values. Therefore, A reaction is faster than others if its propensity is of several orders of magnitude larger than the other propensity values (see the list of abbreviations and notations at the end).

Chemical master equation
In this paper, we consider a network of biochemical reactions at a constant volume.
The network consists ofÑ ≥1 different species fS 1 ; …SÑg. They are spatially homogeneous and interact through M ≥1 reaction channels in thermal equilibrium. The number of counts of each different species defines the state of the system. If all the species are bounded by S, then the approximate number of states in the system would be SÑ [13]. Each state X ≡ ðx 1 ; …xÑÞ T . xÑ, denotes the number of molecules (counts) of each species. For every state, X, the probability satisfies the following CME [8], where P (t) (X) = the probability function, representing the time-evolution of the system, given that t ≥ t 0 and the initial probability is, P ðt 0 Þ ðX 0 Þ, M = elementary chemical reaction channels R 1 , … .. R M , a μ = chemical reaction propensity of channel μ = {1, 2, …. M}, and v μ = the stoichiometric vector that represents a change in the molecular population of the chemical species due to the occurrence of one R M reaction. The system transitions to a new state: X + v μ records the changes in the number of counts of different species when the reactions occur.
We note that a μ (X − v μ )dt is the probability for state (X − v μ ) to transition to state X through chemical reaction, R M , during [t, t + dt), and P M μ¼1 a μ ðXÞdt is the probability for the system to shift from state X as a result of any reaction during dt. If X J ¼ fX 1 ; ……:X SÑ g is the ordered set of possible states of the system indexed by {1, 2…K} having SÑ elements, then Eq. (2) represents the set of ordinary differential equations (ODEs) that determines the changes in probability density P (t) = (P (t) (X 1 ), … P ðtÞ ðX SÑ ÞÞ T .
Once X J is selected, the matrix-vector form of Eq. (2) is described by an ODE: where the transition rate matrix is A = [a i,j ]. If each reaction leads to a different state, X i 0 , then the elements in submatrix A i,j are given as: This equation represents the infinitesimal generator of the Markov process [10,14,15]. Rows and columns are ordered in lowercase letters, i and j respectively. The entry of a i,j of the matrix determines the propensity for the chemical system to transition from one state to another state, given that i ≠ j, are non-negative. The diagonal terms of the matrix are defined by a jj , when i = j and the matrix has a zero-column sum, so its probability is conserved. From Eq. (3) we can derive the P ðt f Þ probability vector at the final time, t f , of interest given an initial density of P ðt 0 Þ : where the matrix exponential function is defined by the convergent Taylor series as [16,17] exp t f A However, algorithms, such as in [13,[18][19][20] truncate Eq. (6) infinite summation to approximate Eq. (3) at the cost of a truncation error.

Initial value problem
If v μ or v M , for μ or M = {1, 2, …. M} are the stoichiometric vectors for R M reaction channels, then we will define the stoichiometric matrix for the system by V μ or V M = [v 1 ; v 2 ; ……v μ ] T . If φ is the sample space and X 0 ∈ φ is the initial state of the system, X J denotes the only set of states in φ. To solve P (t) (X) in Eq. (2) for X ∈ φ, we define the P (t) vector as(P (t) (X)) X ∈ φ or ðP ðtÞ ðXÞÞ X∈X J for a finite set of states, then ∂P ðtÞ ∂t is defined as a vector ð ∂P ðtÞ ∂t Þ X∈φ . Solving the CME involves finding the solution of the initial value problem over a time period using the differential equation Eq. (3) when t > 0, whereas, P ðt 0 Þ is the initial distribution at t = 0. Here, the sample space φ can be infinite for large biochemical systems. Finding the solution for Eq. (3) for the given parameters with a finite set of states X J is a major problem for CME's because in large biochemical systems the size of A will be extremely large.
As seen in Eq. (5), solving Eq. (2) becomes a problem when the model's dimensions grow due to the increase of species present in the system. This is particularly true for large biochemical models. The approximate estimate of SÑ shows how the size of the problem increases. This explosion in size is known as the curse of dimensionality [9,13]. The CME solution given in Eq. (5) has two major parts: (a) the expansion of the state-space, and (b) the approximation of the series. For the expansion of state-space, Finite State Projection (FSP) [21] and Sliding Windows (SW) [18] are used to find the domain. Methods like Krylov subspace [13] and Runge Kutta [22] are commonly used for approximation (of the series) of the CME Eq. (5).
Although CME has been employed and solved explicitly for relatively small biological systems [13, 18-20, 23, 24], computationally complaisant but accurate solutions are still unknown for most significant systems and for large systems which have an infinite (or a very large) number of states. This lack of closed-form solutions has driven the system biology research towards Monte-carlo Algorithms (MC) [25] to capture dynamics. One algorithm, the Stochastic Simulation Algorithm (SSA) by Gillespie [9], has been used in the CME. Although the original FSP state-space expansion has been used in research [21,26], it has some drawbacks [21]. The FSP [21] and its variants [20,24,26,27] are based on r-step reachability [26]. While SW [18] is also a FSP based method, it employs a stochastic simulation algorithm (SSA) to find the domain. This is more effective than FSP and suitable for stiff problems. Add-on weighting functions like GORDE [28] and likelihoods [24] methods are used to improve the expansion. FSP GORDE [28] removes the states with small probabilities before the calculation of Eq. (5). This practice saves computational time, meaning that FSP GORDE performs faster than conventional FSP r-step reachability. However, removing the probabilities before the calculation of Eq. (5) increases the steps error and affects the accuracy of the final solution at t f regardless of whether the state-space is small or large. If one is interested in solving stiff and/or large systems, it will greatly affect the solution.
The FSP variant, Optimal Finite State Projection (OFSP), [20] based on r-step reachability, performs better in terms of producing optimal order domain. It is faster than both FSP and FSP GORDE. However, it is infeasible to use SW for large CME problems because creating hyper-rectangles is a very difficult task. At least four-times the number of SSA simulations are required to minimize the error by half, because of very low convergence rates of routines in MC . The original SSA takes a long time, because one simulation may have several different R M . Recently, the SSA's efficiency has been greatly enhanced by researchers through various schemes such as τ leaps (adaptive) [29,30]. Thus, we compare the OFSP and SSA (τ leaps adaptive) against the ISP in terms of finding the domain, accuracy and computational efficiency. Key to solving the CME remains in finding the right projection size (domain) for large models which would then ensure efficient approximation.
In this paper, we focus primarily on developing the expansion strategy, namely the Intelligent State Projection (ISP) method, to mitigate several problems: the accuracy of the solution, the performance of the method and projection size. The ISP has two variants: the Latitudinal Search (LAS) and the Longitudinal-Latitudinal Search (LOLAS). It treats the Markov chain of a biochemical system as a Markov chain tree structure and states as objects of class node. Based on the dimension of the system, search is performed in a latitudinal way for different model sizes using the ISP LAS method. Whereas, bidirectional search is applied using ISP LOLAS, which quickly expands the state-space up to a specified bound limit. To support the expansion strategy, we also develop the Bayesian Likelihood Node Projection (BLNP) function, based on Bayes' theorem [31,32]. It is adjoined with the ISP variants to determine the likelihood of events at any interval at the molecular population level. BLNP provides confidence to the expansion strategy by assigning probability values to the occurrence of future reactions and prioritizing the direction of expansion. The ISP embedding BLNP function inductively expands the multiple states with the likelihood of occurrence of fast and slow reactions [12]. It also defines the complexity of the system by predicting the pattern of state-space updation, and the depth of the end state from the initial state. When used for any size of biological networks, LAS' memory usage is proportional to the entire width of expansion; it is less than ISP LOLAS. Both methods are feasible and differentiated for various types of biological networks. However, the computational time for both variants depend on the nature of the model and the size of the time step used. At any point, the amount of memory in use is directly proportional to the neighboring states reachable through a single R M reaction. ISP LOLAS uses considerably less memory, even when it retracts to the initial node to track new reactions, then revisiting the depth many times.

Results
Having discussed the CME solution, we now discuss the modeling and integration of the biochemical reaction systems for the ISP methods, as well as the assumptions underlying these methods. Using ISP, we tested its ability to reproduce the model to measure dynamics of the key parameters in the models. The ISP method is a novel, easy-to-use, technique for modeling and expanding the state-space of biochemical systems. It features several improvements in modeling and computational efficiency.
The computational experiment (initializing and solving the model) was conducted on the carbon-neutral platform of Amazon® Web Service Elastic Computing (EC2), instance type large (m5a), running on HVM (hardware virtual environment) virtualization with variable ECUs. We used multicore environment 16vCPU @ 2.2GHz, AMD EPYC 7571 running Ubuntu 16.04.1 with relevant dependencies, and 64GB memory with 8GB Elastic Block Storage (EBS) type General Purpose SSD (GP2) formatted with Elastic File System (EFS). The performance mode was set to General Purpose with input-output per second (IOPS = 100/3000). We used the type bursting throughput mode (see Supplementary Information (SI) 1).

Intelligent state projection
The main aim of the proposed algorithm is to expand the X K iteratively, such that X K contains a minimum number of states carrying the maximum probability mass of the system. To create the sample space for ISP, a Markov chain tree Ѭ [33] was used to visualize a biochemical system to exhibit the transition matrix as directed trees [10,11] of its associated graph. Additionally, the Markov chain tree Ѭ generates sample space of the system to represent Markov processes associated with the Markov chain and the transition matrices of biochemical reaction networks. In following section, we visualize the Markov chain of the biochemical system as a Markov chain graph (tree) for ISP compatibility.

Markov chain as a Markov chain tree
We define the Markov chain tree, Ѭ, [33] as infinite and locally finite. It is a special type of graph with a prominent vertex called a parent node without loops or cycles. If graph G mc is a state-space of the finite state Markov chain with fPðX i ; X i 0 Þ j X i ; X i 0 ∈ G mc g transition probabilities meeting the condition Markov chain tree is a combination of valued G mc random variables with the distributions inductively defined from PðX i ; X i 0 Þ with an initial state, X i ∈ G mc . That being the case, it is easy to expand this class of Markov random field through a Markov chain tree structure for biochemical systems. Furthermore the Markov chain tree and the Markov processes can be equated as explained in [34] for the stochastic analysis.
Since we are interested in aperiodic states in the expansion of state-space, we shall assume the reducibility or simplification of the G mc ; namely for each X i ; X i 0 ∈G mc through Ѭ. Therefore, let us concentrate on the case where G mc is considered as a locally finite connected graph. The transition probabilities of each state are not equal due to the propensities and parameters of different reactions in the biochemical system. Consequently, a Markov chain tree, Ѭ, can be used to visualize a biochemical system process to exhibit a transition matrix as directed trees of its associated graph [10,11]. It can also be used to generate a sample space for the system to represent the Markov processes and the transition matrices of biochemical reaction networks. We discuss the details needed to represent Markov models on trees and working with graphs for state-space later.
If X J is the finite set of cardinality {1, 2……. K} of a Markov chain Ѫ c , then A is the transition probability matrix associated with X J . A state-space is, substantially, a class of a set of states containing the unique state of the system. The arcs between the states represent the transitions from the initial state to the end state. This transition is defined as transient and communicating class in graphs. When all the transitions are combined, every state-space takes the form of a graph and creates the state-space of the system, as shown in Fig. 1

below.
We can now associate chain Ѫ c with the directed graph G mc = (X J , V μ ), where V μ = [v 1 ; v 2 ; ……v μ ]. v μ defines the transition from state X i to X i 0 and is denoted as v μ ¼ fðX i ; X i 0 Þ ; a i; j > 0g. For every transition ðX i ; X i 0 Þ∈X J , then weight ωðX i ; X i 0 Þ is a i, j .
Suppose G mc has a cycle, which starts and terminates at some state, X i ∈ X J . If there is a transition from X i to X i 0 , we add a unique transition by creating a cycle from X i back to itself and then consider the original transition from X i to X i 0 . This contradicts the uniqueness of the walk in tree [35]. In terms of the CTMC of a biochemical system process, the change in molecular population is defined by a stoichiometric vector, so, in G mc , there must be at least one intermediate state that will send the system back to the previous state to create the cycle. This process categorizes the forward and backward reactions given the initial state, X 0 , of the system. The transient class of the transition leads the system to a unique state that defines the forward reaction in the system. In contrast, the communicating class of a transition defines the reversible reaction in the system. We define such systems as transient class systems and communicating class systems. Large biochemical systems are usually a combination of both classes.
A biochemical system is visualized as a tree Ѭ [33] to enable the expansion of the state-space. A tree, Ѭ, is a special form of graph in data structure constituting a set of nodes and a collection of edges (or arcs), each of which connects to an ordered pair of nodes. G mc is considered a directed tree, Ѭ. It is rooted with N 0 = (X 0 , ƌ l ) if it contains a unique walk to N i = (X i , ƌ l + 1) and does not contain any cycles. Meanwhile, if X i ∈ X J \{X 0 } has exactly one outgoing transition away from X 0 it is called an arborescence. If it makes its transition towards N 0 = (X 0 , ƌ l ) it is called an anti-arborescence. An arborescence is a subset ⊆V μ that has one edge out of every node that contains no cycles and has maximum cardinality. For example, if set U = {5, 7, 8, 10} contains 4 elements, then the cardinality of |U| is 4.
If X i and X i 0 are the states other than the initial X 0 state, there is a transition from X i to X i 0 , so X i has at least one transition. Now, suppose X i has two walks, (X i , X i + 1 ) and (X i , X i + 2 ). Concatenating these walks to the walks (X i + 1 , X i 0 ) and (X i + 2 , X i 0 ), respectively, we have two distinct changes in state from X i to X i 0 in G mc . However, in Ѭ, this concatenation is not considered, which makes them Directed Acyclic Graphs (DAG) (see Fig. 2). Most of the biochemical models G mc can be visualized as DAGs irrespective of the nature of the reactions present in the model. Figure 2 shows the equivalent G mc tree of shown in Fig. 1. The trees are less complex as they have no cycles, no self-loops. Yet they are still connected, meaning they can depict the state-space.
The weight of the tree containing all e edges is defined by , where ωðeÞ ¼ ωðX i ; X i 0 Þ ¼ a i; j is the weight of an edge starting from X i and ending at X i 0 when [36]. For systems which have both forward and backward reactions, if n J is the total number of nodes indexed by {1, 2…K} the same as states, then n K is the set of nodes carrying X K , and n 0 K is the set of nodes carrying X 0 K given N 0 root node of the tree Ѭ, then the walk from one node to another node is given by: Ѭ is formed by superimposing the forward transitions between the states X i and X i 0 , with the reverse orientation. Where X i 0 and X i , indicate backward reactions, these are graphically denoted as an individual edge from N i = (X i , ƌ l ) to N i' = (X i , ƌ l + 11) to N i = (X i , ƌ l + 2) in a tree. The N i of ƌ l + 2 can be renamed as a new node. N i + 1 , remains as it is at a different depth from the N i of ƌ l but contains the same state X i . In the expansion, repeated states are not considered in the domain; therefore, any node which carries a similar state is considered the same, regardless of the level and indexing. Consideration of trees for the state-space expansion in ISP not only helps to reduce the complexity but also improves the accuracy of the solution of Eq. (5) by identifying nodes which carry probable states. If the Markov chain graph starts in state X i ∈ X J , then the mean number of transits to any state X i 0 converging to a X i ;X i 0 is given by the (i, i ′ ) th value of U is the set of all arborescences. Let U X i ;X i 0 is the set of all arborescences which have a transition from X i to X i 0 and kU X i ;X i 0 k is the sum of the weights of the arborescences in U X i ;X i 0 then according to the Markov chain tree theorem [33], a X i ;X i 0 is probabilistic in nature. This nature is not only restricted to the systems which have irreducible Markov chains, in which graph G mc is strongly connected while carrying probable state-spaces, but also for the systems that can be simplified by converting to a Markov chain tree and then by reducing that tree by ignoring the states which have low probabilities in space φ.

Expansion criterion for state space
As previously mentioned, the states are indexed using {1, 2……. K} in the domain denoted by set X J . To derive the time, based on the state-space expansion conditions, the probability exponential form of the CME Eq. (5) is evaluated for approximation up to the desired final time t f in steps. To focus on the probable states that contribute most to the probability mass in the domain, we first define the set of non-probable states (those which have the least probability mass) as X 0 K , which are to be bunked. The number of states will usually be infinite, without selecting probable states for the domain. By doing this we can avoid recalculating the probabilities and decrease the computational efforts by keeping the domain small. This bunking can also be applied to the  Fig. 1. This is a special form of graph which has no cycle and no self-loops. It depicts the state-space of the system in the form of a tree (DAG) initial distribution of the system at t 0 . If submatrix A 0 j contains the non-probable set X 0 K of states, then the probability of set will be, The criterion for defining the non-probable states is determined by the τ m tolerance value. A 0 j will only be considered to have non-probable states if, Similarly, submatrix A j has a probable set X K of states if P (t) (X K ) ≥ τ m otherwise, the For any iteration, if P ðtÞ ðX 0 K Þ ≥ τ m then (from Eq. (11)) some states from X 0 K return to X K in the next iteration to increase the accuracy of the approximate solution (Ӕ). The column sum of the approximate solution (Ӕ) of these states is defined as: where, I = (1, …1) T is of an appropriate length. Declaring some states as nonprobable and removing them before calculating the probabilities as seen in [28] will decrease the accuracy of Ӕ with the cumulative step errors. This can be validated from the state probabilities that have been ignored in the domain: We define the step error in terms of the probabilities bunked. If e rror ∝P ðtÞ ðX 0 K Þ then, Every expansion step explores at least one new state and change {X K } but not neces- is satisfied. For ideal systems with a given initial probability of P ðt 0 Þ ðX 0 Þ , the fX 0 K g should be null and so P ðt f Þ ðX 0 K Þ ¼ 0. For such systems fX K g; fX 0 K g ϵ fX J g for final projection and, (18) is the solution of Eq. (3) after the state-space is expanded to X K . However, for large biochemical systems, Eq. (18) may not hold completely true, due to the nature ((fast (R M(fs) ) and slow (R M(sr) )) of some reactions present in the system; therefore, the condition in Eq. (11) will pass the states from X 0 K to X K . The states with the lowest probabilities will be bunked when: This improves the solution. Removing without calculating the probabilities of some states is one of the lags of the current methods [18,20,21,24,[26][27][28]; it is a result of achieving the truncated domain and saving computation time. To address this, we set a P ðtÞ ðX 0 K Þ leakage point based on: where, τ m (leak) for systems will reform Eq. (16) as: which would then zip the X 0 K further by leaking the highest probabilities to X K so that the probability sum is no longer conserved. The motivation of setting this ration is to reconsider (up to 40% of X 0 K ) the bunked states to improve the Ӕ solution and decrease the expansion step error. While modeling the biochemical system, if the slow and fast reaction [12] criterion is considered during expansion, then τ m (leak) will be defined as, We consider Eq. (21) criterion for all the computational experiments in this study. The conditions in Eqs. (21) and (22) will lead to an optimal set of states as, at t d in the domain. When X K is updated at every t step before reaching t f , this creates several intermediate domains which we define as Bound. At t 0 , the domain only has the initial state of the system; therefore, we define the Bound as: ð24Þ After a single t step of expansion, if X K is updated with a new state or set of states, it creates: Here, ƌ l denotes the depth level of the latest state or set of states that has been added in the domain to form Bound upper . This Bound upper is re-labeled and considered as Bound lower for the next t step of the expansion. If the expansion is to be limited in the number of Bounds, then every count(ƃ limit ) leads to: where, ƃ limit is the bound limit. For example, if ƃ limit = 2, then the count(ƃ limit ) will be from 0 → to 1 → to 2. If the count(ƃ limit ) is increased to ƃ limit for I tr th iterations, then Bound upper in the current iteration will be Bound lower for the next iteration. Every Bound lower state will be the strict subset of every consecutive Bound upper given as: and the upper bound as: where Z is the number of Bounds (or intermediate domains). The 2D pyramid domain in Fig. 3 graphically illustrates the increase in the population of states in the domain with the increase in iterations (I tr ). The apex of the pyramid represents the initial state X 0 of the system at Bound lower (1) at t 0 , whereas the base represents the deepest level where the system ends with the final domain carrying set X K with the maximum probability mass.
For large biochemical systems, the number of created Bounds are based on I tr . They have million/billions of states. Expansion can be terminated by defining time t f at which the solution is required. To have an auto break-off point in the expansion, it is first necessary to define the criteria that limits I tr or when no more new states can be searched. Therefore, we define this criterion in the following section. This criterion also applies to biochemical systems which have fast and slow reactions [12].

Cease of criterion after updating
In every expansion step, the domain is validated by Eq. (21) and new states are added in X K as long as: is satisfied for probable states and stops if Eq. (29) is not satisfied. This leads to a point at t f where e rror < τ m , but the expansion can be extended to meet accuracy requirements by re-considering the criteria as: before steps to t f . However, the size of X K obtained through Eqs. (30), (31) and (32) at t f will be greater compared to the size of X K obtained by Eq. (29) at t f , as the latter will have fewer states. In Eqs. (30), (31) and (32), with the increase in the size of A j , the value of the left-hand side will also increase, resulting in an improvement in Ӕ. When considering any Markov process of a biochemical system of any size in which the probability density expands according to Eq.

Computational experimental results
The ISP method is initialized and parameterized using the initial conditions of the models. Due to a large number of mathematical operations and equations, simultaneous parameter predictions with a limited number of experimental values is often complicated for dynamic systems. Therefore, the consistency with the available experimental data was ensured at each step of the ISP. This method has led to the successful development of several functions that integrate large number of processes supporting extensive expansion of the state-space.
To demonstrate the ISP LAS algorithm, we first consider the catalytic reaction system [37] defined by the reactions depicted as a network in Fig. 4 as: In this biochemical system (dimension = 5), reactant P will transform into product E via complex B when reactant S acts as a catalyst for the reaction and produces C. We rewrite this catalytic reaction system as a network of three reactions: with the initial copy counts S 0 = 50, P 0 = 80, B 0 = C 0 = E 0 The reaction rate parameters are k 1 = 1, k 2 = 1000, k 3 = 100. These species counts are used as a state-space to define the model and these copy counts are tracked as ð½S; ½B; ½C; ½P; ½EÞ∈Ñ≔ðx 0 ; x 1 ; x 2 ; x 3 ; x 4 Þ.
In reaction R 1 , the copy count of S is reduced by 1, which increases the copy count of B by 1. In reaction R 2 the copy count of B is reduced by 1, which increases the copy count of C by 1. In contrast, reaction R 3 decreases the B and P counts by 1 and increases the B and E counts by 1. As in R 3 , B acts as a catalyst to convert P to E and B is retained in the same reaction. We can now define the transitions associated with R 1 , R 2 , R 3 in the stoichiometric vector V M matrix as: For LAS method compatibility, the associated Markov chain of this model is converted into a Markov chain tree with the states in terms of the nodes with additional information such as the number of R M reactions required to reach the state. In the growing Markov chain tree, the transition between the nodes: is defined in the typical form of the dictionary Dict. We express the propensity functions of the three reactions in terms of the states ð½S; ½B; ½C; ½P; ½EÞ∈Ñ . Node N 1 = (X 0 , ƌ 1 ) carries the initial state X 0 of the system at an initial depth of level 1. Further, n J = (X K , ƌ 1, 2, … ) is expanded and the states updated by following the LAS order. The corresponding propensities Δa i,j are updated in the A i,j matrix in every iteration, based on the LAS updating trend (for example, see SI 2). The system began with S 0 = 50, P 0 = 80. Gradually all the reactants are transformed to products, E and B. The system ends in n J = (X 1, 2, ……14666 , ƌ l ).   Table 1. The system's conditional probabilities based on its species are shown in Fig. 6.
In three test runs, the ISP LAS run time for the catalytic system was 4677 secs when solving Eq. (5) with 14666 states. The probability of the species in Figure SI 17 (see SI 9) shows the nature of the reactions affecting each species' count in the system. The involvement of species B in all the reactions results in its highest probability at t f . Species B also acts as a catalyst for R 3 , converting species P to E; therefore, both have equal probabilities at the time of solution. Figure 6 shows the total probability bunked at t ′ while progressing with the expansion. Bunking produces an error (w.r.t approximation), with time when the number of states increases with the expansion. LAS produces minimal error of order 10 −5 , as given in Table 1.
To demonstrate the ISP LOLAS algorithm, we consider the coupled enzymatic reactions defined by the reactions depicted as a network in Fig. 7 as: This biochemical system (dimension = 6) describes two sets of enzymatic reactions transforming species S into species P and transforming species P back into S. We rewrite C reactions system as a network of six reactions: with initial copy counts S= 50, E 1 = 20, E 2 = 10, C 1 = C 2 = P= 0 and reaction rate parameters of k 1 = k 4 = 4, k 2 = k 5 = 5, k 3 = k 6 = 1. These species counts are used as a state-space to define the model. These copy counts are tracked as: As in the previous example, we can now define the transitions associated with R 1 , R 2 , R 3 , R 4 , R 5 , R 6 in the stoichiometric vector V M matrix as: For the LOLAS method, the associated Markov chain of this model is converted to a Markov chain tree with the states in terms of nodes with additional information, such  6 Total probability of states bunked at t ′ from the domain of the catalytic system produced by ISP LAS iteration while expanding and solving the CME as the number of R M reactions required to reach the state. In growing Markov chain tree, the transition between the nodes: is defined in the typical form of the dictionary Dict. We express the propensity functions of the six reactions in terms of the states ð½S; ½E 1 ; ½C 1 ; ½P; ½E 2 ; ½C 2 Þ∈Ñ. Node N 1 = (X 0 , ƌ 1 ) carries the initial state X 0 of the system at the initial depth (level 1). Then n J = (X K , ƌ 1, 2, …… ) is further expanded and the states updated by following the LOLAS order. The corresponding propensities Δa i,j are updated in the A i,j matrix in every iteration, based on the given LOLAS updation trend (for example, see SI 2). Initially, the system started with S= 50, E 1 = 20, E 2 = 10. Gradually all reactant species were transformed into products resulting in the system ending in n J = (X 1, 2, ……8296 , ƌ l ).  Table 2. The system's conditional probabilities based on species are shown in Figure SI 18 (see SI 10) Fig. 7 Coupled enzymatic reactions network. The figure shows sixÑ¼6 species, S, E 1 , C 1 , P, E 2 , C 2 , in a network, defining reactions, as given in Eqs. (39) and (40) In three test runs, ISP LOLAS ′ run time for the dual enzymatic reaction network was ≈1614 secs when solving Eq. (5) SI 10) shows the nature of the reactions which affect each species' count in the system. At t f , the probabilities of E 2 and C 2 remain high compared to E 1 and C 1 at different molecular counts. This results in a low probability of P compared to S. We know that this network transforms species S into species P and then transforms species P back into S. Based on the current probabilities of the species at t f , the future probability of P will increase. S will remain the same or decrease. With this change, the probabilities of E 2 and C 2 decrease in comparison to E 1 and C 1 . Figure 9 shows the total probability bunked at t ′ while progressing with the expansion. The bunking produces an error (w.r.t approximation) with time when the number of states increases with the expansion and provided that, LOLAS produces a minimal error of order, 10 −5 , as given in Table 2.
We extend the application of our ISP method to simulate a large model of the G1/S network [38] under the condition of DNA-damage. We want to determine the number of states at different points in time and predict the conditional probabilities of the protein species based on events leading to the formation of different complexes in the system.
The G1/S model (dimension = 28) with a DNA-damage signal transduction pathway is considered to be very stiff in nature, so while some molecular counts of certain proteins increase very rapidly others do so slowly. This makes it tough to solve, even for a short time period. The model is solved for t f = 1.5 sec with ƃ limit = 1, τ m = 1e − 6, t step = 0.1. The systematic exploration of nodes carrying probable states are undertaken in a similar way as discussed in Table SI 4  The nodes are expanded up to t f to enable identification of the reaction channels responsible for variations in the proteins. From the transitioning factor of the 2 nd -tier, we  can see that every node has an average of at least ≈97 possible child nodes carrying states. Further, Dict is expanded for n-tiers of the child nodes to add more states to the domain. Additionally, n J = (X K , ƌ l = 1,2,…… ) is expanded and updated, as per the ISP LOLAS trend (see Table SI 5  The corresponding propensities, Δa i, j , are updated in the A i, j matrix in every iteration, based on the ISP LOLAS update trend. The system started with the initial state of the protein species and gradually, as protein levels change in the system, it exploits the copy counts that shift the system to a new state. The change in protein levels causes the system to transform into new states: here we see the manifestation of the Markov process of the system. The ISP LOLAS captures this process and defines several  Table 3. The size of bound created in each duration reveals that for every step, the growth of the domain is eight-to-ten times the previous size of the domain. The set of nodes N 1 , N 2 , . . …N 3409900 carries unique states representing the set of state(n 3409900 ) = (X 0, 1, 2, …3409899 ) that forms the state-space of the model. It is important to note that some proteins are synthesized and promoted by the network itself, as evidenced by some reactions of the pathway, which increase the frequency of the repeated states. However, ISP LOLAS validation does not consider them for the domain. Equation (5)'s solution, up to t f for the domain, created by ISP LOLAS, is shown in Table 4.
Over three test runs, the ISP LOLAS′ run time for the G1/S model was 1372 secs for solving Eq. (5), with the optimal domain having 3409899 states. The ISP LOLAS response given in Fig. 11, shows the system's probabilities bunked at t ′ during the expansion (w.r.t approximation), when the number of states increases with the expansion, and provided that ISP LOLAS produces minimal errors of the order of 10 −6 , as given in Table 4 and Fig (a) of Fig. 11. We set the checkpoint to examine the initial state's probability over time. The response in Fig (a) of Fig. 11 indicates that the probability of the system remaining in the initial (normal) state decreases with time in the presence of DNA damage, which triggers the change in protein levels.
The conditional probabilities of the species' systems are given in Fig. 14 Figure SI 19 (see SI 11). This leads to the progression to the S-phase, followed by the temporary suspension of cell cycle progression. The increase in probability of x 24 (p53) shows cell support to repair the DNA damage. The parameters and the probabilities relating to x 14 (p21) and x 24 (p53) become important in the case of DNA  Fig. 11 The ISP LOLAS′ response for total probability bunked at t ′ from the domain and checkpoint for examining the initial state probability over time. (a) shows how ISP maintains accuracy by keeping low errors. (b) shows the decline in the probability of the system remaining in the initial state in the presence of DNA damage damage. When combined, the conditional probability of these parameters indicates the involvement of the DNA-damage signal in the transition of G1/S.

Discussion
In this section, we discuss ISP performance, focusing specifically on the speed and accuracy of the expansion, domain size and accuracy of the solution in comparison with other methods.

Comparison with other methods
An approximation of 10 −5 is used to find the approximate number of realizations required by the SSA for the 10 −4 global error. Realizations were computed until the difference was less than 10 −4 between the known distribution and the empirical distribution. Approximately 10 6 and 10 5 runs were required to obtain the right distribution for the catalytic and dual enzymatic reaction networks, respectively. In the catalytic system, we observe (see Table 5) that both versions of ISP are faster than the OFSP of r-step reachability and the SSA of sliding windows. We attribute this greater efficiency to LOLAS having fewer states and less computational time than the OFSP method. LOLAS has better accuracy at t f . Similary, the ISP was much faster than the SSA, and the total number of realizations required from the SSA to have an error at t f still large than that of LOLAS is 10 6 . In the dual enzymatic network, we observe (see Table 5) that both versions of ISP are faster than the OFSP of r-step reachability and the SSA of sliding windows; we attribute the improvement to both ISP variants having an efficient domain with a small approximation error and less computational time than that of the OFSP method and better accuracy at t f . Similarly, both ISP variants were much faster than SSA, as the total number of realizations required to have an empirical distribution with the error at t f is ≈12 times more than the domain produced by ISP.
We also compared the error at t f to determine the solution's efficiency. As seen in the results, the increase in the step error in OFSP affected the solution at t f . Figure 12 (see Fig (a) and (b)), compares the ISP (LAS and LOLAS) with OFSP on the basis of the approximation error at t during the expansion of the catalytic and dual enzymatic reaction networks, respectively. Addressing the step error in ISP and the selection of the probable states results in a more efficient solution at t f compared to OFSP.
The typical firing nature of reactions in the catalytic system makes them stiff. Therefore, the selection of states becomes difficult to approximate. While some species in the system increase abruptly, others do so very slowly because the kinetic parameters (k 1 = 1, k 2 = 1000, k 3 = 100) have large differences: this triggers reactions at different rates. Reaction R 1 , is categorized as a slow reaction in the network: it affects the fast reaction, R 2 . As the computation results of Table 5 show, the ISP found that only 13089 probable states were required to solve the system up to t f . This not only saves computational time (see Fig. 13) compared to OFSP and SSA, and improves the solution's accuracy. In OFSP, applying the compression at every step or after a few steps is still computationally expensive for a model like the catalytic reaction system, as seen in Table 5 and Fig  (a) of Fig. 13.
The network shown in Fig. 7 consists of two interlinked enzymatic reaction systems. These systems transform species S and P into each other via the other species, making the system stiff in nature. The selection of states thus becomes difficult for approximation. This is due to some species (S and E 1 ) in the system increasing abruptly, while others take longer to increase. Some of the kinetic parameters (k 1 = k 4 = 4, k 2 = k 5 = 5) have large differences from other kinetic parameters (k 3 = k 6 = 1): this triggers reactions at different rates. Categorized as the fastest reaction in the network, R affects species S, C, E. It is followed by other reactions involving other species. As the computation   Table 5 show, the ISP LAS indicated that only 8282 probable states are required to solve the system up to t f . Likewise, ISP LOLAS identified that only 8296 probable states are required to solve the system up to t f This saves computational time (see Fig. 13), compared to OFSP and SSA, and improves the solution's accuracy. In OFSP, applying the compression at a defined step or after a few steps is still computationally expensive for models like the dual enzymatic reaction system, as seen in Table 5 and Fig (b) of Fig. 13.
The total computation effort required at every step, when compressing the number of states up to t f , is approximately equal to the total computation effort required when the compression is applied in the gaps in some steps on a set of states up to t f . Moreover, the state-space will remain the same at t f , regardless of when the compression is applied.
A comparison of the computational times in Table 5 shows that both versions of ISP are significantly faster than other methods. Figure 14 shows the CPU utilization (%) of LOLAS and OFSP with respect to run-time (minutes). The dedicated throughput (see SI 1.1) between EC2 and EBS was not used to solve the model. The average CPU exertion is about 60%, which is a considerable workload for a given model. The expansion and approximation began when CPU use was at ≈1.6422% in the catalytic reaction system and ≈1.23% in the dual enzymatic reaction system, at t= 0 sec. It increases up to 60.0% and then drops down to zero at t f .

Theoretical interpretation of methods
Although, SSA recognizes the support and wastes no time in searching for the right domain and creating independent realizations which can be run parallel on multi-core environment, solving the system via Eq. (5) is quicker than creating realizations via stochastic simulation [39][40][41]. This is because the N-term approximation [42,43] of the probability distribution to create the required number of realizations is always less than, or equal to, the minimal support approximation up to same error. These realizations were computed until the difference was less than the prescribed approximation between the known distribution and the empirical distribution. In terms of the system dimension, which is usually defined by the number of species in the system, the approximation of Eq. (5) in ISP becomes smaller problem to solve compared to approximation through SSA. This enables ISP to perform better.
In contrast, OFSP creates a hyper-rectangle and applies truncation to guarantee the minimal order domain for the approximation. OFSP truncates the state space after every few steps to ensure the minimum size of the domain and enable greater computational speed. However, differences in reaction firing changes the probability of some states at a later stage; therefore, truncating the state space in OFSP after every few steps or at every step would remove probable states from the domain, which in turn would affect the accuracy of the solution. As a result of this, OFSP ′ s overall performance is compromised. In contrast, ISP first explores the states based on guided exploration through the BLNP function (see method section (a)) and then leaks the set of states X 0 K which have the lowest probabilities in the bunker at t ′ without removing them (see Eq. (20)). It recalls these sets of states when the probabilities of these states increase at later time.
In ISP, the time and space complexity (refer to SI 7) of removing and accessing the states in the bunking and recalling process is optimum, compared to the overall time and space complexity of the truncation step in OFSP [20]. As seen in Table 5, the number of states present in the domain for the catalytic reaction network in ISP LAS is approximately equal to number of states present in the domain produced by OFSP. Additionally, the number of probable states in the domain for the dual enzymatic reaction system in ISP LAS is quite more as compared to the domain produced by OFSP. However, better complexity and the guided selection of probable states for the domain produces low approximation errors and means that ISP LAS performs better overall than OFS. Similarly, ISP LOLAS outperforms OFSP in finding the optimum domain due to its bi-directional exploratory nature (see methods section (c)). This feature helps ISP LOLAS to achieve a more accurate solution (see Table 5 and Fig. 12) as well as a quicker computational time (see Fig. 13). These benefits are also due to fact that ISP visualizes the state-space as a Markov chain graph or a tree (see Markov chain as a Markov chain tree section) which ultimately decreases the complexity in the expansion phase.

Conclusions
This paper has introduced a novel approach, ISP,to model biochemical systems. This new approach addresses both performance and accuracy problems in CME solutions. Provided all probable states are not added into the domain, up to the desired t f , variants of ISP (LAS and LOLAS) provide systematic ways of expanding the state-space. We have demonstrated the effectiveness of our methods with several experiments using real biological models: the catalytic reaction system, the dual enzymatic reaction system, and the G1/S model (large model). The results and the algorithm's responses reveal improvements in how different sized biological networks can be modeled: even statespaces with 3409900 nodes (see Table 3) carrying states up to ≈3.5 million can be explored within a reasonable time. The results also show that the domain laid out by ISP had an optimal order and was successful in finding probable system states, all the while maintaining high levels of accuracy and efficient computational timing.
We have compared the ISP results against two popular methods: OFSP (r-step reachability) and SSA (τ leaps adaptive). ISP outperformed the other methods, in computational expense, accuracy and projection size. The ISP was more effective in terms of predicting the behavior of the state-space of the system and in performance management, which is a vital step in modeling large biochemical systems. Unlike other methods, the ISP keeps the lowest states probabilities in the bunker without removing (as removed in OFSP) them, before calculation (as removed without calculation in FSP GORDE). It computes the probabilities at t without computing large numbers of realizations (as done in SSA).
The diverging nature of the ISP response, with respect to OFSP in Fig. 19, also shows that the solution improved with t and at a higher t f . For example, in the large model (case study 2), the computation time was 1372 sec and the solution was 3.52e − 06 at t f , which was lower than the small model results (the catalytic reaction system). These results show ISP's compatibility with the distinct size of biochemical models.
These examples have demonstrated that ISP is a very promising technique for system's biology. For stiff models, such as the G1/S and Candida albicans models, the ISP yielded plenty of information. Likewise, it provided opportunities for stochastic analysis of large models. ISP can be used to compute the probabilities of the species up to the required time. One could also use ISP to conduct robustness and sensitivity analysis on the dynamics of biochemical systems and to keep track of what reactions are more active in the system at a particular time. ISP is also able to determine the complexity of the system by defining the bounds with number states and keep track of the nested state-space patterns (called the ISP model blueprint) that were updated at the end of each step. Outlining the patterns of expansion of states to predict the projection folds can be used for updating the new states.
We anticipate that the current structure of the ISP variants can be employed for different classes and varieties of biological systems. Additionally, they can be used to compute the configurations with many reactions, as long as the notable part of the statespace density is present between Bound(Z) lower and Bound(Z) upper . When there was a high probability of the molecular population of the species undergoing several excursions in a fraction of time, then the ISP uses a small t step to capture these moments. While such computations were still challenging in the expansion phase for typical models they can be addressed more closely in combination with the second part of the CME solution: that is, the approximation phase. There are several methods which can be used to address these challenges.
Approximation methods, such as the Krylov sub-space, can be used to effectively compute the matrix exponential times of a vector. While it was mathematically attractive to aggregate the states or decompose the large sparse matrix into a small dense matrix using the Krylov sub-space, this method may not be computationally efficient in the absence of an efficient domain. Performance can also be enhanced by employing the fast math functions, compatible with the multicore environment. We have clearly outlined the core ideas behind the ISP variants. We have highlighted the differences and similarities between them and other methods that cover the computational and theoretical considerations that are essential before any of the approximation methods becomes feasible for an efficient CME solution.

Methods
To understand and predict the dynamics of the state-space response in biochemical systems, we have developed an analytical numerical method called ISP. This method integrates the reactions' propensities describing the Markovian processes through set of nodes governing set of states of the system. The two variants of ISP, named LAS and LOLAS, consist of several modules that incorporate sets of inputs and functions within several compartments. Figure 15 depicts all the ISP modules. The integrated form is discussed later in Tables 7 and 8.
These modules and sub-modules constitute the ISP method. They track key changes in the components that follow changes in the reaction propensities by population and activation of the species. The modules also describe the dynamics of the biochemical system. The method also permits the time form quantification of state-space based on the size and model dimension.
The ISP states expansion strategy is based on the Artificial Intelligence (AI) standards [44][45][46][47], state-space search and relative successor (S uc ) operator or function which perform operations on inputs. AI refers to the study of intelligent agents [48] of a system that perceives and takes action to successfully achieve goals. Most of the problems can be formulated as searches. They can be solved by reducing to one of searching a graph or a tree. To use this approach, we must specify the successor operator which defines the sequence of actions leading from initial to goal state at different time intervals, that lead to the solution.
In terms of AI, we define the state-space as a set of states in the system we can get to by applying S uc to explore new states of a biochemical network. S uc can be applied explicitly, which maps the current state to the output states, or defined implicitly, in that it acts on the current state and transforms it into a new state. In the state-space graph for biochemical networks, we do not define the goal state (or end state) explicitly. Instead, it is defined by S uc implicitly in intervals based on the nature (fast, slow, reversible and irreversible) of the reactions in the system, the duration of expansion and the introduction of stochasticity into the system. This should systematically expand the state-space from X K at t to X K + 1 at t + 1 by going through each node n J at depth ƌ l of the Markov chain graph to evaluate the Markov processes; the expansion aims to occupy most of the probability mass during [X K + X K + 1 ], and the Markov processes can be solved for probability distribution at t + 1.
Where X J is the finite set of states and G mc = (X J , V μ ) is the Markov chain graph on X J associated with A = [a i,j ], given X 0 as the initial state and X K as the set of the explored state, where X 0 ∈ X K then the implicit successor is defined as, Equation (44) gives the new states of the system, where, V μ is the set of stoichiometric vectors v μ function defining the state transitions from any present state X i ∈ X K to new state X i 0 ∉X K . The sample space in the graph contains the unique state of the system stored in a transition matrix, which satisfies Eq. (7) conditions. This transition matrix is a compressed row format (CSR) [49,50] based on an index of rows ⟶ columns delimited by commas generating the dictionary Dict of the model which defines the transitions between nodes in the state-space and the mapping of states. Through S uc , we can know nothing more than the neighbors (child nodes) of the current node (states reachable through a single reaction). We then consider these neighbors (child nodes) as our only goal states; there can be many in numbers. In a situation such as this, search trails are referred to as blind or uninformed searches. In the following section, we discuss the infrastructure of an uninformed search, the type of data structure we will be dealing with.

Infrastructure for searching
A data structure is required to retain the search track in the graph for problem statespace expansion. For each node, N i , of the tree, we create a structure consisting of five elements: (1) N i .State: represents state X i in the state-space corresponding to N i ; (2) N i .Parent: represents the parent node of the child node N i ; (3) N i .Depth: represents the depth of state state X i ; (4) N i .Cost: represents the cost Ͼ N i ;N 0 i of the transition from N i to N 0 i in the statespace; (5) N i .Action: represents the action applied via S uc on the parent node to reach N i .
To explore new states in the system, we consider the initial state state(N 1 ) = (X 0 , ƌ l ) as input to the successor, S uc . Once the expansion is initiated, the Dict will temporarily (in run-time) store the information for the transition from one node to another in the state-space that binds to the reaction propensities a μ . This shift is denoted by an arrow →, which shows multiple transitions from the parent nodes to the child nodes containing the end state. The set of nodes n J ¼ fN 1 ; N 2 ; …::N SÑ g is a data structure that incorporates the Markov chain graph G mc . We explore all the nodes that store the set of states X K as well as some additional information about the state, such as the depth and transition cost, from one state to another in the system. If a set of states(n J ) = X J , then Ͼ N i ;N 0 i is the transition cost to reach stateðN 0 i Þ ¼ X i 0 from state(N 1 ) = X 1 and depth(n J ) = ƌ l defines the depth of the set of nodes in G mc , then the standard relation between a set of nodes and a set of states is given by n J = (X J , ƌ l ) or and the standard relation between a single node and a single state is given by if the transition cost is considered.
For example, Fig (a) of Fig. 16 shows the Markov chain graph, G mc , with n J = 10, ƌ l = 4. Its equivalent tree Ѭ is shown in Fig (b) of Fig. 16 with n J = 15, ƌ l = 5. In the tree nodes N 1 = N 11 = N 12 carries the same state, X 1 at ƌ l = 1, 2 and 3, respectively, where walk N 2 → N 11 and N 7 → N 12 represent the backward reaction of the forward reaction represented by walk N 1 → N 2 and N 1 → N 7 , respectively.
The set of nodes with states are represented as ð45Þ ð46Þ In general, the transition cost, , is defined as: is the summation of all the propensities a μ of the R M reactions that take the system to its final state. For example, to expand to state(N 10 ) = X 10 of Fig (a) of

Fig. 16 is given by
If these are the possible paths for the expansion that expands X K at every iteration then will be defined by the only path that has the lowest P ðtÞ ðX 0 K Þ. This can be generalized as follows: which means that in order to minimize the expansion cost for the optimal domain X K at least one path should have states with high probabilities for X K . It is best to follow the path with , which leaks the minimum probabilities of the system.
For large biochemical models, there exist infinite cases when the node is unreachable from the initial or another node; such cases are ignored when because some probabilities are always dropped in the approximation. Therefore, as defined by the lowest P ðtÞ ðX 0 K Þ is strictly limited to, Upon expanding the root node N 1 , we expand the child nodes carrying new states, and then the child-child nodes are explored. The walk between nodes N i → V μ ðX K ðtÞÞ N i + 1 is defined by dictionary Dict. This represents the occurrence of R M reactions through M elementary channels. For Fig (a) of Fig. 16, the typical form of dictionary is given below: and is indexed with the propensities, [a i,j ], for all the R M reactions. As the propensities are changing by Δa i,j , we consider the recent values of a i,j in every iteration of ISP that corresponds to the reactions involved. To make the feasible for any type of biochemical system (stiff, non-stiff) to capture probable states, it is important to consider the expansion cost for small t step (time step). This may be because there are some cases when to reach two or more different child nodes are equal or very close to each other. In addition, we intend to expand the state-space in the direction of carrying states with high probability mass. To achieve this, we treat or convert our uninformed search to an informed search infrastructure at run-time to have intuitive knowledge beyond our reach. Figure 17 shows the limits of our visibility in the state-space. Consequently, it is important to track the reactions which have high propensity function values. As it is difficult to determine the direction of the expansion, in the following section (a), we develop the post successor function on Bayes' theorem [31,32] to prioritize the expansion direction based only on those reactions that can be triggered at

Bayesian likelihood node projection function
Bayesian methods [32,51] are based on the principle of linking prior probability with posterior probability through Bayes' theorem [31,32]. Posterior probability is the improved form of prior probability, via the likelihood of finding factual support for a valid fundamental hypothesis. Therefore, we employ the standards of Bayes' theorem to develop a function targeted to ensure the quality of the expansion based on R M reactions active in the network at any particular moment. For a concise definition for the purpose of fundamentals, refer to SI 6.
To improve the quality of expansion through a projection function, one may find it useful to remove the set of states which have low probabilities before calculating Eq. (3). However, removing these states will compromise accuracy as the step error will increase at every t. Moreover, removing these probabilities will greatly affect the solution, as defined at t f (at which a solution is required), for large dimension systems which have large state-spaces, as the step error will be much higher due to dropping probabilities without solving Eq. (3). In large systems, any species may change its behavior after a certain number of firing of reactions triggering inactive reactions in the network that will affect the probabilities of the states. If a change in behavior increases the probabilities of certain states, then removing them in an earlier stage is not wise.
Through the Bayesian Likelihood Node Projection (BLNP) function, we seek to predict the posterior probability based on the parent state's probability and calculate the likelihood of the occurrence of reactions that will take the system from the present state to the future state. Through BLNP, we can capture knowledge about the system that will help us to make better predictions about the future state. We are also able to ensure the accuracy of the solution and an optimal domain. Fig. 17 Limits of our visibility in the state-space before expansion. Visualized using a Markov chain graph, where is the initial node and are nodes that are directly reachable from the initial node when exactly one R M occurs. When a further R M occur, the system jumps to other nodes It is important to decide on the direction of the expansion when choosing the future state of the system, as any reaction can occur and take the system to any new state. To understand this situation more clearly on a node level, we assume a Markov chain graph as shown in Fig. 18 of this system which has almost the same number of species count. In Fig. 19, the expansion is at intermediate position, as the initial state state(N 0 ) = X 0 is already expanded and now the expansion of state(N 2 ) = X 2 can be undertaken. To calculate the likelihood of the occurrence of reactions R 1 , R 2 , R 5 , we consider the propensities a i, j as a parameter. Δa i,j depends on the kinetic parameter of the reaction. To assign weight to our belief, we deduce a function that will calculate the probability of reactions occurring and prioritize the expansion in order from reactions resulting in states with high probabilities to reaction giving states with low probabilities. It is important to note that none of the probabilities will be removed before the calculation of Eq. (5). With this function, the likelihood of occurrence of R M can be computed.
We consider each node as a junction of the prior reactions fR To calculate the likelihood of the reactions, it is necessary to have prior information about the occurrence of reactions. If the expansion is to be done at the initial node state(N 0 ) = X 0 (at level 1), then the prior likelihood value b 0 N;N 0 is considered as the initial probability or as ≈1. Once the initial node has been explored, we can calculate the likelihood of the reactions inductively. To calculate the probabilities b 1, N′ , …. , b N, N′ of the occurrence of R 1 , …. , R M , we first calculate the weighted probabilities P N, 1 (ω), … … .,P N;N 0 ðωÞ of a system leaving any state by: where, Once bðN N1;::NM jb 0 1;N…:N 0 ;N Þ is calculated for all the adjacent nodes, the values are arranged in descending order. Every value is bound to one reaction and represents the likelihood of the occurrence of the reaction that takes the system from the present node to the child nodes. Based on likelihood values (highest to lowest), the corresponding reactions are considered one by one and labelled as true events for expansion. For example, if a system has R 1 , R 2 , R 3 reactions that bound to BLNP likelihood values in order from highest to lowest, respectively, then three events take the system to new state. When R 1 is considered for expansion, R 2 and R 3 are labeled false events and R 1 as the true event. When the second highest BLNP likelihood value is considered, which is for R 2 , then it is labeled the true event and the others, R 1 , R 3, are labeled false events. Similarly, the last and lowest BLNP likelihood value is for R 3 , which is labeled as the true event and the others as false events. All states are added in the domain in order from the 1 st true event to the  Figure 18 shows the Markov chain tree for selection present at level 2 (assuming that the initial node is already expanded). Here we calculate the weighted probability of a system leaving state(N 2 ) = X 2 by: similarly, P 2, 4 (ω) = 0.3333 and P 2, 5 (ω) = 0.3418. At level 2, the conditional probability of the occurrence of reaction R 1 , given the probability of occurrence of reaction R 1 at level 1, is given by: Similarly, the occurrence of reaction R 1 at level 2, given the probability of occurrence of reaction R 6 at level 1, is given by: If at level 1, state(N 1 ) = X 1 and at level 2, state(N 2 ) = X 2 are explored through R 1 then we say that this is a true event and temporarily consider other events false events with respect to the other reactions. Such a condition holds true for the other two cases, when, at level 1, state(N 1 ) = X 1 is explored through R 1 followed by an exploration of state(N 2 ) = X 2 either by R 2 or R 5 . Given b  Table 6. The likelihood values of future reactions cannot be equal, as they are based on the probabilities of occurrence of prior reactions.
From Fig. 18 and Table 6, we can infer, based on the prior reactions for R M , where M = {1, 6} that: The BLNP function cannot be used standalone for expansion because it only assigns weightage to direction for expansion. In the Intelligent state projection section, we have derived the condition for our expansion strategies to work with the Markov chain graph state-space and defined the criteria for the formation of bounds (domain formed at anytime t) with time. The BLNP function (with expansion strategies), will choose the probable states in large biochemical systems where it is important to capture the moments at time t that define a system's behavior. BLNP will be useful for identifying the most active reactions in the system while guiding the expansion towards the set of states with high probability mass.

Latitudinal search strategy
We delve deeper into the first subroutine of the ISP called the Intelligent State Projection Latitudinal Search (ISP LAS). Figure 20 manifests the infrastructure of the LAS strategy, showing G mc , the queue and the domain. LAS' queue data structure is based on the FIFO (First In, First Out) method. In this method, the oldest state added to the queue is considered first. We define and exploit the direction of expansion step-by-step based on intuitive knowledge (as discussed in section (a)), gained from the probability of future reaction events. We follow the conditions (as discussed in the results section). At level ƌ l the states are expanded only after all the states at level ƌ l − 1 have been expanded; that is, the search is undertaken level-by-level and depth ƌ l increases in every I tr iteration. In the case of networks with reversible reactions, the ISP condition will prevent LAS from returning to the state it came from and also prevent transitions containing cycles resulting in DAG with no repetition of any state whatsoever; however, changes in propensities a i,j are validated. Verifying the explored states in X K in iterations ensures that the algorithm completes and that a deadlock in the state transition cycles cannot occur.
The time complexity of LAS depends on the average transitioning factor Ŧ and depth ƌ l and is given by (see SI 7 for detailed discussion), For the nodes at the deepest level ƌ l , all walks are valid except for the very last node which stores the end state of the system. Therefore, once the end state is found, based on Eq. (20), LAS will zip X 0 K , further leaking the highest probabilities to X K for the solution of Eq. (3) which includes the end state of the system. As no state is ever repeated in the domain, space complexity will decrease when the set of states X 0 K is bunked at t ′ seconds in iterations if Eqs. (19) and (20) are satisfied. In Eq. (13), P ðtÞ ðX 0 K Þ is computed according to Eq. (5) (the exponential form of the CME), where τ m is the tolerance and I is the identity matrix. Due to this stepping bunking of X 0 K from X K , the time complexity O(Ŧ d + 1 ) reduces to , where |X J | is the size of the state-space [13]. In contrast, the expansion of new nodes carrying similar states tend to increase however, repetitive states are ignored.
If the input τ m is too small, the algorithm automatically uses the default value of sqrt(eps). Here sqrt is the square root and eps is the default value of the epsilon on machine. The expansion of child nodes containing state(N i ) = X i stops if the condition of Eq. (32) is not satisfied. If the criterion of slow and fast reaction [12] is considered, then the condition of Eq. (31) or (32) is used, depending on the number of R M(sr) and R M(fs) . Table 7 shows the steps of the LAS method with the embedded BLNP function, from steps 4a to 5b.
LAS will be optimal if the transitions between all the states are uniform; that is, all the R M reactions have the same propensity values. However, in real biochemical models, this condition is unusual. To see a step-by-step demonstration of the ISP LAS algorithm on a toy model, refer to SI 2. We now turn our attention to the second variant of the ISP. We apply the method to a toy model to see how it differs from LAS.

Longitudinal-latitudinal search strategy
Here, we delve deeper into the second sub-routine of ISP called the Intelligent State Projection Longitudinal Latitudinal Search (ISP LOLAS). Figure 21 visually represents the infrastructure of the LOLAS strategy, showing the G mc , stack and the domain. The stack data structure of LOLAS is based on the LIFO (Last In, First Out) method. In this method, the newest state added to the stack is considered first. In particular, we define the bound limit and exploit the direction of the expansion step-by-step based on intuitive knowledge (as discussed in section (a)), gained from the probability of future reaction events and follow the conditions (as discussed in the Results section). Furthermore, we show step-bystep how nodes carrying states are explored in a bidirectional way and how these states were updated in the domain in I tr iterations.
The states at level ƌ l are expanded only after the neighboring states at level ƌ l − 1 have been expanded for R M : that is, the search is undertaken level-by-level. Depth ƌ l increases Table 7 Steps of ISP latitudinal search (LAS) algorithm Step 0: Inputs: Initial node N 0 , a μ , v μ , tol τ m , t f , t step Initialize: Step 1: Start from parent node N i = (X 0 , ƌ l ) ← Current State of the system at t d , Step 2: Flag the current node as explored, update A and add the state X i in the domain so that; if 1 − I T exp (t. A j ). P (t) (X 0 ) ≥ τ m (leak) holds true go to Step 3; else stop the algorithm Step 3: Sort exp(t. A j ). P (t) (X 0 ) and shift the set of states in X 0 K at t 0 having smallest probabilities, if Step 4a: Extend the graph dictionary Dict by v μ (X i (t)) by 1 level to check all the nodes n j = (X j , ƌ l , Ͼ Ni ;N 0 i ð minÞÞ adjacent to N i : Bound upper ← R M (Bound lower ) reachable by exactly R M reactions (from fast to slow) having Ͼ Ni ;N 0 i ð minÞ. If n K = (X K , ƌ l , Ͼ Ni ;N 0 i ð minÞÞ be the set of adjacent nodes such that n K ∈ n J then go to next Step, Step 4b: Compute the BLNP function for n K ∈ Bound upper : bðN N1;::NM jb  in the same I tr iteration up to a certain ƃ limit (bound limit). The expansion limit is set by ƃ limit , ƌ step (depth step). In contrast to LAS, it is not set by depth ƌ l ,. The LOLAS search updates the dictionary Dict of G mc by the stoichiometric vector function, v μ (X(t)) on state at level ƌ l to explore the child nodes carrying states on levels ƌ l + 1, ƌ l + 2 … .. ƌ l + l, where l = {1, 2…. ∞} and then retracts to level ƌ l at which new state exploration decisions can be made. In the case of networks with reversible reactions, the ISP conditions will prevent LOLAS from returning to the state it came from and prevent transitions containing cycles resulting in DAG with no repetition of any state whatsoever; however, the change in propensities a i, j is validated. Verifying the explored states in X K in iterations ensures that the algorithm completes and that deadlocks in the state transition cycles cannot occur.
In the absence of ƃ limit , the algorithm will not retract. It will explore longitudinally by tracking only one R M reaction. In addition, the algorithm will not terminate with an optimal order domain carrying a maximum probability mass. This would lead to an increase in the approximation error. Instead, it will terminate when carrying only those set of states as a result of tracking only a few R M , creating an insufficient domain for approximation. Therefore, by default, the value of ƃ limit ≥ 1 is kept for large systems and can be increased depending upon the model's dimension and the availability of the testing environment's random access memory (RAM). LOLAS' worst-case time complexity depends on the average transitioning factor Ŧ. Depth ƌ l is given by (see SI 7 for a detailed discussion): ð57Þ Table 8 Steps of ISP longitudinal latitudinal search (LOLAS) algorithm Step 0: Inputs: Initial node N 0 , ƌ step , ƃ limit , a μ , v μ , tol τ m , t f , t step Initialize: Step 1: Initialize count(ƃ limit ) and start from parent node N i = (X 0 , ƌ l ) ← Current state of the system at t d , Step 2: Flag the current node as explored, update A and add the state X i in the domain so that; if 1 − I T exp (t. A j ). P (t) (X 0 ) ≥ τ m (leak) holds true go to Step 3; else stop the algorithm.
Step 3: Sort exp(t. A j ). P (t) (X 0 ) and shift the set of states in X 0 K at t 0 having smallest probabilities, if P (t) (X K ) ≥ τ m (leak) > P ðtÞ ðX 0 K Þ and at t d update X K ←X K − X 0 K Step 4a: For ƌ step , extend the graph dictionary Dict by v μ (X i (t)) for count(ƃ limit ) to check all the nodes n j = (X j , ƌ l , Ͼ Ni ;N 0 i ð minÞÞ adjacent to N i : Bound upper ← R M (Bound lower ) reachable by exactly R M reactions (from fast to slow) having Ͼ Ni ;N 0 i ð minÞ. If n K = (X K , ƌ l , Ͼ Ni ;N 0 i ð minÞÞ be the set of adjacent nodes such that n K ∈ n J , then go to the next Step, Step 4b: Compute the BLNP function for n K ∈ Bound upper : bðN N1;::NM jb 0 1;N:…N Step 5a: If n K = (X K , ƌ l , Ͼ Ni ;N 0 i ð minÞÞ∈domain, then update the values of the set of states X K present in domain and take domain ← domain previous ∪ domain and go back to Step 1; else If n K = (X K , ƌ l , Ͼ Ni ;N 0 i ð minÞÞ∉domain, then add it to the stack in order, according to reachability and go to next Step, Step 5b: Step 6: Pop of the top nodes n K = (X K , ƌ l , Ͼ Ni ;N 0 i ð minÞÞ from the stack and add the set of states X K in the domain as domain ← domain + X K and take domain previous ∪ domain, and go to next Step, Step 7: If count(ƃ limit ) = ƃ limit creates Bound upper = {domain} up to ƃ limit then label Bound lower ← Bound upper and go back to Step 1; else if count(ƃ limit ) < ƃ limit creates {domain} up to count(ƃ limit ) then go to next Step, Step 8: count(ƃ limit )← count(ƃ limit ) + 1 and go to Step 4a Output: domain with probable states LOLAS only stores the transition path to the end state besides the neighbors of each relevant node in the exploration. Once all descendants are updated with the relevant propensities in the projection, it discards the node from the domain (explored), making it ready for the approximation. LOLAS first considers the R 1 reaction and the corresponding stoichiometric vector v 1 of the system, to explore all the neighboring states up to bound limit ƃ limit . It then considers R 2 , R 3 , … ..,R M for the same ƃ limit and the corresponding v 2 , v 3 , … .. , v M to explore the states. For count(ƃ limit ), LOLAS retracts to the R 1 reaction and explores the new neighboring states longitudinally. It then reconsiders R 2 , R 3 , … .. , R M to explore the other states in a similar fashion. Provided with this reaction tracking pattern, the BLNP function alters this trend and guides this tracking by considering reactions in a different order based on their propensities and the number of probable states of the system.
If the system is ending in a set of state X K carried by n K at t f , then LOLAS will explore the states efficiently, as long as count(ƃ limit ) ≤ ƃ limit , otherwise count(ƃ limit ) is reset for further expansion. Choosing the appropriate ƌ limit and ƃ step depends on the type of biochemical reaction network and the computing configuration. Starting with a depth 1 → ƃ limit , LOLAS explores all the states until they return null. It then resets the count(ƃ limit ) and retracts to explore again. In most cases, fewer states are positioned at the lower level. They increase at a higher level when the number of active R M reactions increases, so retracting provides the ability to track all the reactions simultaneously. The nature of the LOLAS expansion means that it is able to find more states at any time t compared to LAS. It is also able to find them at the deepest level of the graph. The states at depth ƌ l are explored once, the states at depth ƌ l − 1 are explored twice, states at depth ƌ l − 2 are explored three times and so on, until it has explored all the system's states. If the input τ m is too small, the algorithm automatically uses the default value of sqrt(eps). Here sqrt is the square root and eps is the default value of the epsilon on machine. The expansion of the child nodes containing state(N i ) = X i stops if the condition of Eq. (32) is not satisfied. If the slow and fast [12] reaction criterion is considered, then either Eq. (31) or (32) conditions are used depending on the number of R M(sr) and R M(fs) . Table 8 shows the LOLAS method, with an embedded BLNP function from steps 4a to 5b.
Refer to SI 3 for the step-by-step demonstration of the ISP LOLAS algorithm, where we assume the same toy model system.

Additional file 1.
Abbreviations a i,j or a μ : Propensity of chemical reaction; Δa i,j : Change in propensity; a 0 1;N 0 , … .,a 0 N;N 0 : Propensities of the prior reactions; a fr : Probability of a jump process from state X i − 1 to X i per unit time; a rv : Probability of a jump process from state X i to X i − 1 per unit time; A or A i,j : Defines the transition between i, j and its weightage; ƃ limit : Exploration bound limit in LOLAS; b