{'fraction_non_alphanumeric': 0.042787171355114756, 'fraction_numerical': 0.016026422954942107, 'mean_word_length': 4.185675360062281, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reviving the Two-State Markov Chain Approach

Probabilistic Boolean networks (PBNs) is a well-established computational framework for modelling biological systems. The steady-state dynamics of PBNs is of crucial importance in the study of such systems. However, for large PBNs, which often arise in systems biology, obtaining the steady-state distribution poses a significant challenge. In this paper, we revive the two-state Markov chain approach to solve this problem. This paper contributes in three aspects. First, we identify a problem of generating biased results with the approach and we propose a few heuristics to avoid such a pitfall. Second, we conduct an extensive experimental comparison of the extended two-state Markov chain approach and another approach based on the Skart method. We analyze the results with machine learning techniques and we show that statistically the two-state Markov chain approach has a better performance. Finally, we demonstrate the potential of the extended two-state Markov chain approach on a case study of a large PBN model of apoptosis in hepatocytes.


INTRODUCTION
C OMPUTATIONAL modelling plays a prominent role in providing a system-level understanding of processes that take place in a living cell. However, it faces significant challenges when modelling realistic biological systems due to the size of the statespace that needs to be considered. Hence, profound understanding of biological processes asks for the development of new methods that would provide means for formal analysis and reasoning about large systems.
We focus on the steady-state dynamics of biological processes modelled as discrete-time Markov chains (DTMCs). This is the case, for example, when the biological system under study is cast into the mathematical/computational framework of probabilistic Boolean networks (PBNs) [1], [2], [3]. In these or other discretetime models, e.g., dynamic Bayesian networks, the real (considered as continuous) time is not modelled. Instead, the evolution of the system is abstracted as a sequence of consecutive events. These coarse-grained models have been successfully applied in many systems biology studies and proved their predictive power [4]. In fact, for the study of large regulatory systems they remain the only reasonable solution. Extrapolating the ordinary differential equations model of a single elementary building block of the network (e.g., a gene) to the whole large system would result in a prohibitively complex model. Moving towards a higher-level description by ignoring the molecular details allows to grasp the system-level behaviour of the network [5]. In consequence, these coarse-grained formalisms are broadly applied in systems biology. One of the key aspects in the analysis of such dynamic systems is the comprehension of their steady-state (long-run) behaviour. For example, cellular phenotypes were hypothesised to be characterised by a collection of sets of states, referred to as attractors, that the system reaches and in which it remains forever in a long run. Another complementary conjecture is that attractors correspond to functional cellular states such as proliferation, apoptosis, or differentiation [6]. These interpretations may cast new light on the understanding of cellular homeostasis and cancer progression [1]. In this work, we focus on the computation of steady-state probabilities which are crucial for the determination of long-run influences and sensitivities. These are measures that quantify the impact of genes on other genes and that enable the identification of elements with highest impact.
So far the huge-state space, which often characterises dynamical models of biological systems, tempers the application of the above mentioned techniques in the analysis of realistic biological systems. In fact, approximations with the use of Markov chain Monte Carlo (MCMC) techniques are the only viable solution to this problem [7]. However, due to the difficulties with the assessment of the convergence rate to the steady-state distribution (see, e.g., [8]), certain care is required when applying these methods in practice. A number of statistical methods exist which allow to empirically determine when the simulation can be stopped, with the Skart method [9] considered to be the stateof-the-art technique. In this study, we focus on another method, i.e., the two-state Markov chain approach. It was introduced in 1992 by Raftery and Lewis, and Shmulevich et al. [7] proposed its application to the analysis of PBNs in 2003. However, to the best of our knowledge, since then it has not been widely applied for the analysis of large PBNs. In this paper, we aim to revive it by demonstrating its usefulness for approximating steady-state probabilities of large PBNs, which often arise in systems biology as computational models of cellular gene regulatory networks.
This paper contributes in three aspects. First, we identify a problem of generating biased results by the original two-state Markov chain approach, due to the size of the initial sample with which the approach needs to start. We propose three heuristics for the approach to avoid some unfortunate initialisations. Secondly, we perform an extensive evaluation and comparison of the extended two-state Markov chain approach with the state-ofthe-art Skart method on a large number of randomly generated PBNs. We apply techniques originating from machine learning to analyse the obtained results. We show that in the context of gene regulatory networks modelled with PBNs, actually the two-state Markov chain approach seems to perform better than the Skart method in most cases in terms of computational cost. Thirdly, we demonstrate the potential of the two-state Markov chain approach on a study of a large, 91-node PBN model of apoptosis in hepatocytes. In particular, we show that the two-state Markov chain approach facilitates the quantification of both the steady-state influences of genes on their target genes and long-run sensitivities in large networks.

Finite discrete-time Markov chains (DTMCs)
Let S be a finite set of states. A (first-order) DTMC is an Svalued stochastic process {X t } t∈N with the property that the next state is independent of the past states given the present state. Formally, P(X t+1 = s t+1 | X t = s t , X t−1 = s t−1 , . . . , X 0 = s 0 ) = P(X t+1 = s t+1 | X t = s t ) for all s t+1 , s t , . . . , s 0 ∈ S. Here, we consider time-homogenous Markov chains, i.e., chains where P(X t+1 = s | X t = s), denoted P s,s , is independent of t for any states s, s ∈ S. The transition matrix P = (P s,s ) s,s ∈S satisfies P s,s 0 and ∑ s ∈S P s,s = 1 for all s ∈ S. We denote by π a probability distribution on S. If π = π P, then π is a stationary distribution of the DTMC (also referred to as the invariant distribution). A path of length n is a sequence s 1 → s 2 → · · · → s n such that P s i ,s i+1 > 0 and s i ∈ S for i ∈ {1, 2, . . . , n}. State q ∈ S is reachable from state p ∈ S if there exists a path such that s 1 = p and s n = q. A DTMC is irreducible if any two states are reachable from each other. The period of a state is defined as the greatest common divisor of the lengths of all paths that start and end in the state. A DTMC is aperiodic if all states in S are of period 1. A finite state DTMC is called ergodic if it is irreducible and aperiodic. By the famous ergodic theorem for DTMCs [10], an ergodic chain has a unique stationary distribution being its limiting distribution (also referred to as the steady-state distribution) given by lim n→∞ π 0 P n , where π 0 is any initial probability distribution on S. In consequence, the limiting distribution for an ergodic chain is independent of the choice of π 0 . The steady-state distribution can be estimated from any initial distribution by iteratively multiplying it by P.
The evolution of a first-order DTMC can be described by a stochastic recurrence sequence X t+1 = φ (X t ,U t+1 ), where {U t } t∈N is an independent sequence of uniformly distributed real random variables over [0, 1] and the transition function φ : S × [0, 1] → S satisfies the property that P(φ (s,U) = s ) = P s,s for any states s, s ∈ S and for any U, a real random variable uniformly distributed over [0, 1]. When S is partially ordered and the transition function φ (·, u) is monotonic, then the chain is said to be monotone ( [11]).

Probabilistic Boolean networks (PBNs)
A PBN G(V, F ) consists of a set of binary-valued nodes (also referred to as genes) V = {v 1 , v 2 , . . . , v n } and a list of sets j ∈ F i is a Boolean function defined with respect to a subset of nodes referred to as parent nodes of v i . There is a probability distribution on each F i ∈ F : c (i) j is the probability of selecting f (i) j ∈ F i as the next predictor for v i and it holds that ∑ l(i) j=1 c (i) j = 1. We denote by v i (t) the value of node v i at time point t ∈ N. The state space of the PBN is S = {0, 1} n and it is of size 2 n . The state of the PBN at time t is determined by s(t) = (v 1 (t), v 2 (t), . . . , v n (t)). The dynamics of the PBN is given by the sequence (s(t)) ∞ t=0 . We consider here independent PBNs where predictor functions for different nodes are selected independently of each other. The transition from s(t) to s(t + 1) is conducted by randomly selecting a predictor function for each node v i from F i and by synchronously updating the node values in accordance with the selected functions. There are N = ∏ n i=1 l(i) different ways in which the predictors can be selected for all n nodes. These combinations are referred to as realisations of the PBN and are represented as n-dimensional function vectors f k = ( f (1) In PBNs with perturbations, a perturbation parameter p ∈ (0, 1) is introduced to sample the perturbation vector γ(t) = (γ 1 (t), γ 2 (t), . . . , γ n (t)), where γ i (t) ∈ {0, 1} and P(γ i (t) = 1) = p for all t and i ∈ {1, 2, . . . , n}. Perturbations provide an alternative way to regulate the dynamics of a PBN: the next state is determined as s(t + 1) = F t (s(t)) if γ(t) = 0 and as s(t + 1) = s(t) ⊕ γ(t) otherwise, where ⊕ is the exclusive or operator for vectors. The perturbations, by the latter update formula, allow the system to move from any state to any other state in one single transition, hence render the underlying Markov chain irreducible and aperiodic. Therefore, the dynamics of a PBN with perturbations can be viewed as an ergodic DTMC [2]. The transition matrix is given by P s, , where 1 is the indicator function and η(s, s ) is the Hamming distance between states s, s ∈ S. According to the ergodic theory, adding perturbations to any PBN assures the long-run dynamics of the resulting PBN is governed by a unique limiting distribution, convergence to which is independent of the choice of the initial state. However, the perturbation probability value should be chosen carefully, not to dilute the behaviour of the original PBN. In this way the 'mathematical trick', although introduces some noise to the original system, allows to significantly simplify the analysis of the steady-state behaviour.
The density of a PBN is measured with its function number and parent nodes number. For a PBN G, its density is defined as where n is the number of nodes in G, N F is the total number of predictor functions in G, and ω(i) is the number of parent nodes for the ith predictor function.
Within the framework of PBNs the concept of influences is defined; it formalizes the impact of parents nodes on a target node and enables its quantification ( [12]). The concept is based on the notion of a partial derivative of a Boolean function f with respect to variable x j (1 ≤ j ≤ n): where ⊕ is addition modulo 2 (exclusive OR) and for l ∈ {0, 1} x ( j,l) = (x 1 , x 2 , . . . , x j−1 , l, x j+1 , . . . , x n ).
The influence of node x j on function f is the expected value of the partial derivative with respect to the probability distribution D(x): Let now F i be the set of predictors for x i with corresponding probabilities c (i) of node x k on the predictor function f (i) j . Then, the influence of node x k on node x i is defined as: The long-term influences are the influences computed when the distribution D(x) is the stead-state distribution of the PBN.
We define and consider in this study two types of long-run sensitivities.
Definition 2.1. The long-run sensitivity with respect to selection probability perturbation is defined as where · l denotes the l-norm, π is the steady-state distribution of the original PBN, p ∈ [0, 1] is the new value for c (i) j , and π[c (i) j = p] is the steady-state probability distribution of the PBN perturbed as follows. The jth selection probability for node x i is replaced withc The remaining selection probabilities of the original PBN are unchanged.
Definition 2.2. The long-run sensitivity with respect to permanent on/off perturbations of a node x i as where π,π[x i ≡ 0], andπ[x i ≡ 1] are the steady-state probability distributions of the original PBN, of the original PBN with all f (i) ∈ F i replaced byf (i) ≡ 0, and all f (i) ∈ F i replaced byf (i) ≡ 1, respectively.
Notice that the definition of long-run sensitivity with respect to permanent on/off perturbations is similar but not equivalent to the definition of long-run sensitivity with respect to 1-gene function perturbation of [12].

The Two-state Markov Chain Approach
The two-state Markov chain approach [13] is a method for estimating the steady-state probability of a subset of states of a DTMC. In this approach the state space of an arbitrary DTMC is split into two disjoint sets, referred to as meta states. One of the meta states, numbered 1, is the subset of interest and the other, numbered 0, is its complement. The steady-state probability of meta state 1, denoted π 1 , can be estimated by performing simulations of the original Markov chain. For this purpose a two-state Markov chain abstraction of the original DTMC is considered. Let {Z t } t 0 be a family of binary random variables, where Z t is the number of the meta state the original Markov chain is in at time t. {Z t } t 0 is a binary (0-1) stochastic process, but in general it is not a Markov chain. However, as argued in [13], a reasonable assumption is that the dependency in {Z t } t 0 falls off rapidly with lag. Therefore, a new process {Z (k) t } t 0 , where Z (k) t = Z 1+(t−1)k , will be approximately a first-order Markov chain for k large enough. A procedure for determining appropriate k is given in [13]. The first-order chain consists of the two meta states with transition probabilities The split of the state space into meta states is marked with dashed ellipses. (b) Projecting the behaviour of the original chain on the two meta states results in a binary (0-1) stochastic process. After potential subsampling, it can be approximated as a first-order, twostate Markov chain with the transition probabilities α and β set appropriately.
α and β between them. See Figure 1 for an illustration of the construction of this abstraction.
The steady-state probability estimateπ 1 is computed from a simulated trajectory of the original DTMC. The key point is to determine the optimal length of the trajectory. Two requirements are imposed. First, the abstraction of the DTMC, i.e., the twostate Markov chain, should converge close to its steady-state distribution π = [π 0 π 1 ]. Formally, t satisfying |P[Z ε for a given ε > 0 and all i, j ∈ {0, 1} needs to be determined. t is the so-called 'burn-in' period and determines the part of the trajectory of the two-state Markov chain that needs to be discarded. Second, the estimateπ 1 is required to satisfy P[π 1 − r π 1 π 1 + r] s, where r is the required precision and s is a specified confidence level. This condition is used to determine the length of the second part of the trajectory used to computeπ 1 , i.e., the sample size. Now, the total required trajectory length of the original DTMC is then given by The functions m and n depend on the transitions probabilities α and β and are given by where Φ −1 is the inverse of the standard normal cumulative distribution function. For the completeness of the presentation, the detailed derivations of the expressions for m and n are given in the Appendices A and B.
Since α and β are unknown, they need to be estimated. This is achieved iteratively in the two-state Markov chain approach of [13]. It starts with sampling an arbitrary initial length trajectory, which is then used for estimating the values of α and β . M and N are calculated based on these estimates. Next, the trajectory is extended to reach the required length, and α and β values are re-estimated. The new estimates are used to re-calculate M and N. This process is iterated until M + N is smaller than the current trajectory length. Finally, the resulting trajectory is used to estimate the steady-state probability of meta state 1. For more details, see [13]. Notice however the small oversights in the formulas for m (absolute value missing in the denominator) and n (the inverse of Φ should be used) therein.

TWO-STATE MARKOV CHAIN APPROACH: THE INITIALISATION PROBLEM
Shmulevich et al. [7] proposed to use the two-state Markov chain approach for analysing the steady-state dynamics of PBNs in 2003. However, since then it has not been widely applied. We aim to revive the two-state Markov chain approach by demonstrating its usefulness for approximating steady-state probabilities of large PBNs. In this section, we first identify an initialisation problem of the original approach due to the size of the initial sample, this particular problem can lead to biased results. We then propose three heuristics to extend the approach for avoiding unfortunate initialisations.
Given good estimates of α and β , the theory of the twostate Markov chain approach presented above guarantees that the obtained value satisfies the imposed precision requirements. However, the method starts with generating a trajectory of the original DTMC of an arbitrarily chosen initial length, i.e., M 0 + N 0 = 1 + (m 0 − 1)k + 1 + (n 0 − 1)k, where m 0 is the 'burnin' period and n 0 is the sample size of the two-state Markov chain abstraction. An unfortunate choice may lead to initial estimates of α and β that are biased and result in the new values of M and N such that M + N is either smaller or not much larger than the initial M 0 +N 0 . In the former case the algorithm stops immediately with the biased values for α, β and, more importantly, with an estimate for the steady-state probability that does not satisfy the precision requirements. The second case may lead to the same problem. As an illustration we considered a two-state Markov chain with α = 24 11873 (0.0020214) and β = 24 25 (0.96). The steadystate probability distribution was [0.997899 0.002101]. With k = 1, ε = 10 −6 , r = 10 −3 , s = 0.95, m 0 = 5, and n 0 = 1, 920 the first estimated values for α and β were 1 1918 (0.0005214) and 1, respectively. This subsequently led to M = 2 and N = 1, 999, resulting in a request for the extension of the trajectory by 76. After the extension, the new estimates for α and β were 1 1997 and 1, respectively. These estimates gave M = 2, N = 1, 920, and the algorithm stopped. The estimated steady-state probability distribution was [0.99950 0.00050], which was outside the prespecified precision interval given by r. Independent 10 4 runs resulted in estimates of the steady-state probabilities that were outside the pre-specified precision interval 10% of times. Given the rather large number of repetitions, it can be concluded that the specified 95% confidence interval was not reached.
The reason for the biased result is the unfortunate initial value for n 0 and the fact that the real value of α is small. In the initialisation phase the value of α is underestimated and n(α, β ) calculated based on the estimated values of α and β is almost the same as n 0 . Hence, subsequent extension of the trajectory does not provide any improvement to the underestimated value of α since the elongation is too short.
To identify and avoid some of such pitfalls, we consider a number of cases and formulate some of the conditions in which the algorithm may fail to achieve the specified precision. To start, let n 0 be the initial size of the sample used for initial estimation of α and β . Neither α nor β is zero. It might be the case that the initial sample size is not big enough to provide non-zero estimates for both α and β . If this is the case, n 0 is doubled and the trajectory is elongated to collect a sample of required size. This is repeated iteratively until non-zero estimates for α and β are obtained. In the continuation we assume that n 0 provides non-zero estimates for both α and β . Then, the smallest possible estimates for both α and β are greater than 1 n 0 . For a moment, let us set an upper bound value for n 0 to be 10 4 . For most cases this boundary value is reasonable. Notice however that this is the case only if the real values of α and β are larger than 10 −4 . In general, the selection of a proper value for n 0 heavily depends on the real values of α and β , which are unknown a priori. From what was stated above, it follows that both first estimates for α and β are greater than 10 −4 . The following cases are possible.
(1) If both α and β are small, e.g., less than 0.1, then we have that 10 −4 < α, β < 0.1 and n(α, β ) > 72, 765 as can be seen by investigating the n(·, ·) function. In this case the sample size is increased more than 7-fold which is reasonable since the two-state Markov chain seems to be bad-mixing by the first estimates of the values for α and β and the algorithm asks for a significant increase of the sample size. We therefore conclude that the bad-mixing case is properly handled by the algorithm.
(2) Both first estimates of α and β are close to 1. If α, β ∈ [0.7, 0.98], the value of n(α, β ) is larger than 19, 000. If both α, β > 0.98, then the size of the sample drops, but in this case the Markov chain is highly well-mixing and short trajectories are expected to provide good estimates.
(3) The situation is somewhat different if one of the parameters is estimated to be small and the other is close to 1 as in the example described above. The extension to the trajectory is too small to significantly change the estimated value of the small parameter and the algorithm halts. Considering the above cases leads us to the observation that the following situation needs to be treated with care: The estimated value for one of the parameters is close to 1 n 0 , the value of the second parameter is close to 1, and n(α, β ) is either smaller or not significantly larger than n 0 .
First approach: pitfall avoidance. To avoid this situation, we determine n 0 which in principle could lead to inaccurate initial estimates of α or β and such that the next sample size given by n(α, β ) would practically not allow for an improvement of the estimates. As stated above, the 'critical' situation may take place when one of the parameters is estimated to be very small, i.e., close to 1 n 0 , and the increase in the sample size is not significant enough to improve the estimate. If the initial estimate is very small, the real value is most probably also small, but the estimate is not accurate. If the value is underestimated to the lowest possible value, i.e., 1 n 0 , on average the improvement can take place only if the sample size is increased at least by n 0 . Therefore, with the trade-off between the accuracy and efficiency of the method in mind, we propose the sample size to be increased at least by n 0 . Then the 'critical' situation condition is n(α, β ) < 2 n 0 . By analysing the function n(·, ·) as described in details in Appendix D, we can determine the values of n 0 that are 'safe', i.e., which do not satisfy the 'critical' condition. We present them in Table 1 for a number of values for r and s.
Second approach: controlled initial estimation of α and β . The formula for n is asymptotically valid provided that the values for α and β are known. However, these values are not known a priori and they need to be estimated. Unfortunately, the original approach does not provide any control over the quality of the initial estimate of the values of these parameters. In certain situation, e.g., as in the case discussed above, the lack of such control mechanism may lead to results with worse statistical confidence level than the specified one given by s. In the discussed example s = 95%, but this value was not reached in the performed experiment. In order to address this problem, we propose to extend the initial phase of the two-state approach algorithm in the following way. The algorithm samples a trajectory of the original DTMC and estimates the values of α and β . We denote the estimates asα andβ , respectively. Next, the algorithm computes the sample size required to reach the s confidence level that the true value of min(α, β ) is within a certain interval. For definiteness, we assume from now on thatα <β , which suggests that min(α, β ) = α. During the execution of the procedure outlined in the following the inequality may be inverted. If this is the case, the algorithm makes corresponding change in the consideration of α and β .
The aim is to have a good estimate for α. Notice that the smallest possible initial value ofα isgreater than 1 n 0 . We refer to 1 n 0 as the resolution of estimation. Given the resolution, one cannot distinguish between values of α in the interval (α − 1 n 0 ,α + 1 n 0 ). In consequence, if α ∈ (α − 1 n 0 ,α + 1 n 0 ), then the estimated valueα should be considered as optimal. Hence, one could use this interval as the one which should contain the real value with specified confidence level. Nevertheless, although the choice of this interval usually leads to very good results, as experimentally verified, the results are obtained at the cost of large samples which make the algorithm stop immediately after the initialisation phase. Consequently, the computational burden is larger than would be required by the original algorithm to reach the desired precision specified by r and s parameters in most cases. In order to reduce this unnecessary overhead, we consider the interval (α −α 2 ,α + α 2 ), which is wider than the previous one wheneverα > 1 n 0 and leads to smaller sample sizes.
The two-state Markov chain consists of two states 0 and 1, i.e., the two meta states of the original DTMC. We set α as the probability of making the transition from state 0 to state 1 (denoted as 0 → 1). The estimateα is computed as the ratio of the number of transitions from state 0 to state 1 to the number of transition from state 0. Let n 0,α be the number of transitions in the sample starting from state 0. Let X i , i = 1, 2, . . . , n 0,α , be a random variable defined as follows: X i is 1 if ith transition from meta-state 0 is 0 → 1 and 0 otherwise.
Notice that state 0 is an accessible atom in the terminology of the theory of Markov chains, i.e., the Markov chain regenerates after entering state 0, and hence the random variables X i , i = 1, 2, . . . , n 0,α , are independent. They are Bernoulli distributed with parameter α. The unbiased estimate of the population variance from the sample, denotedσ 2 , is given byσ 2 =α · (1 −α) · n 0,α n 0,α −1 . Due to independence,σ 2 is also the asymptotic variance and, in consequence, the sample size that provides the specified confidence level for the estimate of the value of α is given by The Markov chain is in state 0 with steady-state probability β α+β . Then, given that the chain reached the steady-state distribution, the expected number of regenerations in a sample of size n is given by n·β α+β . Therefore, the sample size used to estimate the value of α with the specified confidence level s is given by n α = α+β β · n α,s (α, n 0,α ). As the real values of α and β are unknown, the estimated valueŝ α andβ can be used in the above formula. If the computed n α is bigger than the current number of transitions n 0,α , we extend the trajectory to reach n α transitions from 0 to 1 and re-estimate the values for α and β using the extended trajectory. We repeat this process until the computed n α value is smaller than the number of transitions used to estimate α. In this way, good initial estimates for α and β are obtained and the original two-state Markov chain approach using the formula for n(α, β ) is run. Third approach: simple heuristics. When performing the initial estimation of α and β , we require both the count of transitions from state 0 to state 1 and the count of transitions from metastate 1 to state 0 be at least 3. If this condition is not satisfied, we proceed by doubling the length of the trajectory. In this way the problem of reaching the resolution boundary is avoided. Our experiments showed that this simple approach in many cases led to good initial estimates of the α and β probabilities. Discussions. The first approach provides us with safe initial starting points. As can be seen in Table 1, there might however be no safe starting point in certain conditions. Nevertheless, the first approach can be used in the initialisation phase of the other two approaches. The second approach introduces a new iteration process to provide a good estimate of α or β . The third one modifies the two-state Markov chain approach by adding only one extra restriction and therefore is the most simple one. We have verified with experiments that the last two approaches have the potential to make the two-state Markov chain approach meet the predefined precision requirement even in the case of an unlucky initial sample size. As a small example, we show in Table 2 the results for verifying two PBNs each of eight nodes. For each of the PBNs, we compute the steady-state probability for one subset of states using three different approaches: 1) the original twostate Markov chain approach (columns 'Original'), the proposed second approach (columns '2nd') and the proposed third approach (columns '3rd'). The precision and confidence level are set to 0.001 and 0.95 respectively. We repeat the computation for 1000 times and count the percentage of times that the calculated result is within the precision requirement (shown in columns labelled 'Computed confidence level'). As can be seen in Table 2, the original two-state Markov chain approach fails to meet the confidence level requirement while the both proposed approaches can meet the requirement. Due to its simplicity, we use the third approach in the remaining of the paper.

EVALUATION
In this section, we focus on verifying the performance of the two-state Markov chain approach with another related method called the Skart method [9]. We use the tool ASSA-PBN [14], [15] as the platform for this verification. ASSA-PBN is a tool specially designed for steady-state analysis of large PBNs; it includes the two-state Markov chain approach with the simple heuristics presented in Section 3 and the Skart method. For the steady-state analysis of large PBNs, applications of these two  methods necessitate generation of trajectories of significant length. To make this efficient, we applied the alias method [16] to sample the consecutive trajectory state. This enables ASSA-PBN, e.g, to simulate 12, 580 steps within 1s for a 2, 000 nodes PBN, which is hundreds of times faster than the related tool optPBN [17]. In Section 4.1, we briefly describe the Skart method. We present an empirical comparison of the performance of these two methods in Section 4.2.

The Skart Method
We choose the Skart method [9] as a reference for the evaluation of the performance of the two-state Markov chain approach. The Skart method is a successor of ASAP3, WASSP, and SBatch methods, which are all based on the idea of batch means [9]. It is a procedure for on-the-fly statistical analysis of the simulation output, asymptotically generated in accordance with a steady-state distribution. Usually it requires an initial sample of size smaller than other established simulation analysis procedures [9]. Briefly, the algorithm partitions a long simulation trajectory into batches, for each batch computes a mean and constructs an interval estimate using the batch means. Further, the interval estimate is used by Skart to decide whether a steady state distribution is reached or more samples are required. For a more detailed description of this method, see [9].
The Skart method differs in three key points with the two-state Markov chain approach. First, it specifies the initial trajectory length to be at least 1, 280, while for the two-state Markov chain approach this information is not provided. This difference, however, does not change the fact that the two methods can provide the same accuracy guarantee providing that the unlucky choice of the initial trajectory length of the two-state Markov chain approach is fixed as mentioned in the previous section. Second, the Skart method applies the student distribution for skewness adjustment while the two-state approach makes use of the normal distribution for confidence interval calculations. Thirdly, the twostate Markov chain approach does not require to keep track of the simulated trajectories; instead, the statistics (e.g., the α and β as in Figure 1) of the trajectories are enough. the Skart method, however, requires to keep track of the simulated trajectories, which consumes a large memory in the cases of large size trajectories.

Performance evaluation
To compare the performance of the two methods, we randomly generated 882 different PBNs using ASSA-PBN. ASSA-PBN can randomly generate a PBN which satisfies structure requirements given in the form of five input parameters: the node number, the minimum and the maximum number of predictor functions per node, finally the minimum and maximum number of parent nodes for a predictor function. We generated PBNs with node numbers from {15, 30, 80, 100, 150, 200, 300, 400, 500, 1000, 2000}. We assigned the obtained PBNs into three different classes with respect to the density measure D: dense models with density 150-300, sparse models with density around 10, and in-between models with density 50-100. The two-state Markov chain approach and the Skart method were tested on these PBNs with precision r set to the values in {10 −2 , 5 × 10 −3 , 10 −3 , 5 × 10 −4 , 10 −4 , 5 × 10 −5 , 1 × 10 −5 , 5 × 10 −6 , 1 × 10 −6 }. We set ε to 10 −10 for the two-state Markov chain approach and s to 0.95 for both methods.
The experiments were performed on a HPC cluster, with CPU speed ranging between 2.2GHz and 3.07GHz. ASSA-PBN is implemented in Java and the initial and maximum Java virtual machine heap size were set to 503MB and 7.86GB, respectively. We collected 5596 valid (precision being smaller than probability) results with the information on the PBN node number, its density class, the precision value, the estimated steady-state probabilities computed by the two methods, and their CPU time costs. The steady-state probabilities computed by the two methods are comparable in all the cases (data not shown in the paper). For each experimental result i, we compare the time costs of the two methods. Let t TS (i) and t Skart (i) be the time cost for the two-state Markov chain approach and the Skart method, respectively. We say that the two-state Markov chain approach is by k per cent faster than the Skart method if (t Skart (i)−t TS (i)) t Skart (i) k 100 . The definition for the Skart method to be faster than the two-state Markov chain approach is symmetric. In Table 3 we show the percentage of cases in which the two-state approach was by k per cent faster than Skart and vice versa for different k. In general, in nearly 70% of the results, the two-state Markov chain approach was faster than the Skart method. The number of cases the two-state Markov chain approach was faster than the Skart method is also larger than in the opposite case.
Next, we analyse the results with a machine learning technique, i.e., logistic regression, in MATLAB. We use the node number, the precision, and the density class as features. We label each result as 1 if the two-state Markov chain approach is by k per cent faster than the Skart method and as 0 otherwise. We plot the receiver operating characteristic (ROC) curve, which is commonly used to illustrate the performance of a binary classifier against varying discrimination threshold and we give the computed area under the curve (AUC) for different k in Figure 2a. When k 15, the AUC value is over 0.7, which means that the prediction is very good. In another word, for a given PBN and precision requirement, we are able to predict whether the two-state Markov chain approach will be by 15 per cent faster than the Skart method in a very high accuracy rate.
We show in Table 4 (left part) the regression coefficient estimates of the three features. Clearly, the precision plays an important role in the prediction since the absolute value is always the largest. We further analyse how the performance of the two methods change with precision and show in Table 5 the percentage of cases that the two-state Markov chain approach is faster than the Skart method with respect to different precisions. The two-     state Markov chain approach has a larger chance to be faster than the Skart method in all the studied precisions, especially when the precision r is low (e.g., 10 −2 ) or very high (e.g., equal to or less than 10 −5 ). Notably, the chance that the two-state Markov chain approach is faster than the Skart method becomes very large when the precision is very high. This is due to the fact that the Skart method requires a large memory to keep track of the large size trajectory when the precision is high. The CPU performance drops when operating on a large memory; on the other hand, the Skart method may run out of memory. Moreover, we analyse the situation when the Skart method is by k per cent faster than the two-state Markov chain approach. This time it becomes difficult to make an accurate prediction as the largest AUC is only about 0.72 for k = 20 (see Figure 2b). Besides, the coefficient estimates in the right part of Table 4 also vary a lot with k, and precision is not always the dominating factor. The detailed experiment data can be obtained at http://satoss.uni. lu/software/ASSA-PBN/benchmark/benchmark.xlsx.
From the above analysis, we conclude that the two-state Markov chain approach outperforms the Skart method (state-ofthe-art) in analysing large PBNs, especially for computing steadystate probabilities with very high precision.

A BIOLOGICAL CASE STUDY
In [18], a large-scale Boolean network of apoptosis (see Figure 3) in hepatocytes was introduced, where the assigned Boolean interactions for each molecule were derived from literature study.
In [17], the original multi-value Boolean model was cast into the PBN framework: a binary PBN model, so-called 'extended apoptosis model' which comprised 91 nodes (state-space of size 2 91 ) and 102 interactions was constructed. In this extended version the possibility of activation of NF-κB through Caspase 8 (C8*), as described in [17], was included. The model was fitted to steady-state experimental data obtained in response to six different stimulations of the input nodes, see [17] for details.
As can be seen from the wiring of the network, the activation of complex2 (co2) by RIP-deubi can take place in two ways: 1) by a positive feedback loop from activated C8* and P → tBid → Bax → smac → RIP-deubi → co2 → C8*-co2 → C8*, and 2) by the positive signal from UV-B irradiation (input nodes UV(1) or UV(2)) → Bax → smac → RIP-deubi → co2. The former to be active requires the stimulation of the type 2 receptor (T2R). The latter way requires complex1 (co1) to be active, which cannot happen without the stimulation of the TNF receptor-1. Therefore, RIP-deubi can activate co2 only in the condition of co-stimulation by TNF and either UV(1) or UV(2). In consequence, it was suggested in [17] that the interaction of activation of co2 via RIP-deubi is not relevant and could be removed from the model in the context of modelling primary hepatocyte. However, due to the problem with efficient generation of very long trajectories in optPBN toolbox, quantitative analysis was hindered and this hypothesis could not be verified ( [17]). In this work, we take up this challenge and we quantitatively investigate the relevancy of the interaction of activation of co2 via RIP-deubi. We perform an extensive analysis in the context of co-stimulation by TNF and either UV(1) or UV(2): we compute long-term influences of parent nodes on the co2 node and the long-run sensitivities with respect to various perturbations related to specific predictor functions and their selection probabilities. For this purpose we apply the two-state Markov chain approach as implemented in our ASSA-PBN tool [15] to compute the relevant steady-state probabilities for the best-fit models described in [17]. Due to the efficient implementation, the ASSA-PBN tool can easily deal with trajectories of length exceeding 2 × 10 9 for this case study.
We consider 20 distinct parameter sets of [17] that resulted in the best fit of the 'extended apoptosis model' to the steadystate experimental data in six different stimulation conditions. In [17], parameter estimation was performed with steady-state measurements for the nodes apoptosis, C3ap17 or C3ap17 2 depending on the stimulation condition considered, and NF-κB. The optimisation procedure used was Particle Swarm and fit score function considered was the sum of squared errors of prediction (SSE) and the sum was taken over the three nodes in the six stimulation conditions. We took all the optimisation results from the three independent parameter estimation runs of [17], each containing 7500 parameter sets. We sorted them increasingly with respect to the cost function value obtained during optimisation, removed duplicates, and finally took the first 20 best-fit parameter sets.
As mentioned above, we fix the experimental context to costimulation of TNF and either UV(1) or UV (2). We note that originally in [18] UV-B irradiation conditions were imposed via a multi-value input node UV which could take on three values, i.e., 0 (no irradiation), 1 (300 J/m 2 UV-B irradiation), and 2 (600 J/m 2 UV-B irradiation). In the model of [17], UV input node was refined as UV(1) and UV(2) in order to cast the original model into the binary PBN framework. Therefore, we consider in our study two cases: 1) co-stimulation of TNF and UV(1) and 2) co-stimulation of TNF and UV (2). Node co2 has two independent predictor functions: co2 = co1 ∧ FADD or co2 = co1 ∧ FADD ∧ RIPdeubi. The selection probabilities are denoted as c (co2) 1 and c (co2) 2 , respectively. Their values have been optimised in [17].
We start with computing the influences with respect to the steady-state distribution, i.e., the long-term influences on co2 of each of its parent nodes: RIP-deubi, co1, and FADD, in accordance with the definition in Section 2.2. Notice that the computation of the three influences requires several joint steady-state probabilities to be estimated with the two-state Markov chain approach, e.g., (co1=1,FADD=1,RIP-deubi=0) or (co1=1,FADD=0). Each probability determines a specific split of the original Markov chain. For example, in the case of the estimation of the joint steady-state probability for (co1=1,FADD=0), the states of the underlying Markov chain of the apoptosis PBN model in which co1=1 and FADD=0 constitute meta state 1 and all the remaining states form meta state 0. Therefore, the estimation of influences is computationally demanding. The summarised results for the 20 parameter sets are presented for the co-stimulation of TNF and UV(1) or TNF and UV(2) in Table 6. They are consistent across the different parameter sets and clearly indicate that the influence of RIP-deubi on co2 is small compared to the influence of co1 or FADD on co2. However, the influence of RIP-deubi is not negligible.
We take the analysis of the importance of the interaction between RIP-deubi and co2 further and we compute various longrun sensitivities with respect to selection probability perturbation. In particular, we perturb the selection probability c (co2) 2 by ±5%, i.e., we set the new value by multiplying the original value by (1 ± 0.05), and compute in line with Definition 2.1 how the joint steady-state distribution for (apoptosis,C3ap17,NFκB) differs from the non-perturbed one with respect to the l 1 norm, i.e., || · || 1 . We notice that the computation of the full steadystate distribution for the considered PBN model of apoptosis is practically intractable, i.e., it would require the estimation of 2 91 values. Therefore, we restrict the computations to the estimation of eight joint steady-state probabilities for all possible combinations of values for (apoptosis,C3ap17,NFκB), i.e., the experimentally measured nodes. Each estimation is obtained by a separate run of the two-state Markov chain approach with the split into meta states determined by the considered probability as explained above in the case of the computation of long-term influences. To compare the estimated distributions we choose the l 1 norm after [19], where it is used in the computations of similar types of sensitivities for PBNs to these defined in Section 2.2. Notice that the l 1 norm of the difference of two probability distributions on a finite sample space is twice the total variation distance. The latter is a well-established metric for measuring the distance between probability distributions defined as the maximum difference between the probabilities assigned to a single event by the two distributions (see, e.g., [20]). Additionally, we check the difference when c (co2) 2 is set to 0 (and, in consequence, c (co2) 1 is set to 1). The obtained results for the 20 parameter sets in the conditions of co-stimulation of TNF and UV(1) and co-stimulation of TNF and UV(2) are summarised in Table 7. In all these cases, the sensitivities are very small. Therefore, the system turns to be insensitive to small perturbations of the value of c (co2) 2 . Also the complete removal of the second predictor function for co2 does not cause any drastic changes in the joint steady-state distribution for (apoptosis,C3ap17,NF-κB).
Finally, we compute the long-run sensitivity with respect to permanent on/off perturbations of the node RIP-deubi in accordance with Definition 2.2. As before, we consider the joint steadystate distributions for (apoptosis,C3ap17,NF-κB) and we choose the l 1 -norm. The results, given in Table 8, show that in both variants of UV-B irradiation the sensitivities are not negligible and the permanent on/off perturbations of RIP-deubi have impact on the steady-state distribution.
To conclude, all the obtained results indicate that in the context of co-stimulation of TNF and either UV(1) or UV(2) the interaction between RIP-deubi and co2 plays a certain role. Although the elimination of the interaction does not invoke significant changes to the considered joint steady-state distribution, the longterm influence of RIP-deubi on co2 is not negligible and may be important for other nodes in the network.

DISCUSSION AND CONCLUSION
Most current tools for statistical model checking, a simulationbased approach using hypothesis testing to infer whether a stochastic system satisfies a property, are restricted for bounded properties which can be checked on finite executions of the system. Recently, both the Skart method [21] and the perfect 6: Long-term influences of RIP-duebi, co1, and FADD on co2 in the 'extended apoptosis model' in [17] under the co-stimulation of both TNF and UV(1) or UV(2).   simulation algorithm [22] have been explored for statistical model checking of steady state and unbounded until properties. The perfect simulation algorithm for sampling the steady-state of an ergodic DTMC is based on the indigenous idea of the backward coupling scheme [11]. It allows to draw independent samples which are distributed exactly in accordance with the steady-state distribution of a DTMC. However, due to the nature of this method, each state in the state space needs to be considered at each step of the coupling scheme. If a DTMC is monotone, then it is possible to sample from the steady-state distribution by considering the maximal and minimal states only [11]. This was exploited in [22] for model checking large queuing networks. Unfortunately, it is not applicable to PBNs with perturbations. In consequence, the perfect simulation algorithm is only suited for at most medium-size PBNs and large-size PBNs are out of its scope. Thus, we have only compared the performance of the two-state Markov chain approach with the Skart method.
Moreover, in this study we have identified a problem of generating biased results by the original two-state Markov chain approach and have proposed three heuristics to avoid wrong initialisation. Finally, we demonstrated the potential of the twostate Markov chain approach on a study of a large, 91-node PBN model of apoptosis in hepatocytes. The two-state Markov chain approach facilitated the quantitative analysis of the large network and the investigation of a previously formulated hypothesis regarding the relevance of the interaction of activation of co2 via RIP-deubi. In the future, we aim to investigate the usage of the discussed statistical methods for approximate steady-state analysis in a research project on systems biology, where we will apply them to develop new techniques for minimal structural interventions to alter steady-state probabilities for large regulatory networks.

TIONS
Let {Z t } t≥0 be a discrete-time two-state Markov chain as given in Figure 1b. Z t has the value 0 or 1 if the system is in state 0 or state 1 at time n, respectively. The transition probabilities satisfy 0 < α, β < 1 and the transition matrix for this chain has the following form Matrix P has two distinct eigenvalues: 1 and λ = (1 − α − β ). Notice that |λ | < 1. The chain is ergodic and the unique steady-state distribution is π = [π 0 π 1 ] = [ β α+β α α+β ]. Let E π (Z t ) denote the expected value of Z t for any fixed t ≥ 0, with respect to the steady-state distribution π. We have that E π (Z t ) = α α+β . The m-step transition matrix can be written, as can be checked by induction, in the form where λ is the second eigenvalue of P. Suppose we require m to be such that the following condition is satisfied for some ε > 0. If e 0 = [1 0] and e 1 = [0 1], then for j ∈ {0, 1} we have that where T is the transposition operator. For any vector v = [v 1 v 2 . . . v n ] T ∈ R n we use |v| to denote [|v 1 | |v 2 | . . . |v n |] T . Therefore, condition (1) can be rewritten as For j = 0 and j = 1 the above simplifies to respectively. Therefore, it is enough to consider the following two inequalities which, since α, β > 0, can be rewritten as Equivalently, m has to satisfy By the fact that |λ m | = |λ | m this can be expressed as Then, by taking the logarithm to base 10 on both sides 1 , we have that and in consequence, since |λ | < 1 and log |λ | < 0, log (|λ |) .

APPENDIX B DERIVATION OF THE SAMPLE SIZE
By the Law of Large Numbers for irreducible positive recurrent Markov chains Z n → π 1 a. s. with n → ∞, where Z n = 1 n ∑ n t=1 Z t . Now, by a variant of the Central Limit Theorem for nonindependent random variables 2 , for n large, Z n is approximately normally distributed with mean π 1 = α α+β and asymptotic variance σ 2 as = 1 n αβ (2−α−β ) (α+β ) 3 , see Section C for the derivation of the asymptotic variance. Let X be the standardised Z n , i.e., X = Z n − π 1 σ as / √ n .
If follows that X is normally distributed with mean 0 and variance 1, i.e., X ∼ N(0, 1). Now, we require n to be such that the condition P[π 1 − r ≤ Z n ≤ π 1 + r] = s is satisfied for some specified r and s. This condition can be rewritten as P[−r ≤ Z n − π 1 ≤ r] = s, 1. In fact, by the formula for change of base for logarithms, the natural logarithm (ln), the logarithm to base 2 (log 2 ), or a logarithm to any other base could be used to calculate m instead of log. Notice that m does not depend on the choice of the base of the logarithm! 2. Notice that the random variables Z t , Z t+1 which values are consecutive states of a trajectory are correlated and are not independent. and further as P[−r · √ n σ as ≤ Z n − π 1 σ as / √ n ≤ r · √ n σ as ] = s, which is P[−r · √ n σ as ≤ X ≤ r · √ n σ as ] = s.

MULA DERIVATIONS
We start with analysing the minimum values n(·, ·) can attain. The function is considered on the domain D = (0, 1] × (0, 1] and, as mentioned before, the estimated values of α and β are within the range [ 1 n 0 , 1]. Computing the partial derivatives, equating them to zero, and solving for α and β yields α = −β , which has no solution in the considered domain. Hence, the function has neither local minimum nor maximum on D. Let us fix β for a moment and consider n(α, β ) as a function of α. We denote it as n β (α). By differentiating with respect to α, we obtain