On computation of the steady-state probability distribution of probabilistic Boolean networks with gene perturbation (Q433954)

The probabilistic Boolean network (PBN) with gene perturbation framework supposes that at every state of the PBN, each gene may (independently to the others) be activated (with a probability \(p\)) or inhibited (with probability \(1-p\)) due to a external stimuli. It is known that the network evolves following a Markov chain with the additional useful property to be ergodic if \(p>0\) with transition probability \(\tilde A\) given by \[ \tilde A = (1-p)^n A + \tilde P, \] where \(n\) is the number of genes of the network, \(A\) the transition matrix of the non perturbed PBN and \(\tilde P\) is a matrix depending on the perturbation probability \(p\). The purpose of this paper is to study the relative error bound between the steady-state probability distribution \(\pi\) of the unperturbed network (if any) and the steady-state probability distribution \(\tilde{\pi}\) of the perturbed network. On the other hand, using the useful structure of the transition matrix \(\tilde A\), the authors propose three different algorithms to compute efficiently \(\tilde \pi\). The main result on the perturbation bound is Theorem \(3.3.\) which gives a bound for the relative change in the limiting probability distributions \(\pi\) (if any) and \(\tilde \pi\) with respect to the quadratic norm, namely, a bound for \( \| \pi - \tilde \pi \|_{2} / \| \pi \|_2.\) This bound depends on the largest and the smallest singular values of some explicit matrices obtained from the stationary properties : \(A \pi = \pi \), \(\tilde A \tilde \pi = \tilde \pi\) and \(\| \pi \|_1 = \| \tilde \pi\|_1 =1\). From their main result, the authors deduce a more meaningful result by making the bound depending on the perturbation probability \(p\) and not depending on individual states, namely, \[ \frac{ \| \pi - \tilde \pi \|_{2}}{\| \pi \|_2} \leq \begin{cases} \frac{1-(1-2p)^n}{\sqrt{\theta} \sigma_{\text{min}}^{+}(B)}, & \text{if} \quad 0< p \leq 1/2 \\ \frac{2p}{\sqrt{\theta} \sigma_{\text{min}}^{+}(B)}, & \text{if} \quad 1/2<p<1, \end{cases} \tag{1} \] where \(\theta\) belongs to a bounded interval and \(B\) is some explicit matrix and where \(\sigma_{\text{min}}^{+}(B)\) denotes the smallest positive singular value of \(B\). We remark that general perturbation bounds exist in the perturbation theory of stochastic matrices and most of them are of the form \[ \| \pi - \tilde \pi \| \leq \kappa \| \tilde A -A \| \tag{2} \] for some matrix norm \(\| \cdot \|\), where \(\kappa\) is called the condition number; see, e.g [\textit{P. J. Schweitzer}, J. Appl. Probab. 5, 401--413 (1968; Zbl 0196.19803)] for the original work on this area and [\textit{E. Seneta}, Stat. Probab. Lett. 17, No. 2, 163--168 (1993; Zbl 0777.60065)] for a resume of some perturbation analysis results. It has to be remarked that the result (1) shows in particular that if the perturbation probability \(p\) goes to \(0\) then \(\| \pi - \tilde \pi\|_2\) tends toward \(0\). This is expected since when \(p\) is near \(0\) then there is a high probability (tending to \(1\)) than no gene of the network is perturbed at every state of the evolution of the system, so that we have \(A \approx \tilde A\). Thus, the previous remark make a link between the perturbation bounds (2) and (1). The second part of the work gives numerical tools for computing efficiently the steady-state probability distribution \(\tilde \pi\). The authors use the structure and properties of the matrix \(I-\tilde A\) to propose fast algorithms for the computation of the steady-state probability distribution using the recursive stationarity formula: \[ \tilde A \tilde \pi^{k+1} = \tilde \pi^{k}, \qquad k=1,2, \dots, \tag{3} \] which sequence of solutions \((\tilde \pi^{k})_{k \geq 2}\) semi-converge to the solution of \(\tilde A\tilde \pi = \tilde \pi\), for every initial probability distribution \(\tilde \pi^{1}\). Numerical tests are performed in different cases, depending on the choice of the perturbation probability \(p\) and the number of genes \(n\), and compared (in term of the CPU time, residual errors \(\| \tilde A \tilde \pi^{k} - \tilde \pi^{k} \|_2/ \| \tilde A \tilde \pi^{0} - \tilde \pi^{0} \|_2\) and the number of iterations for a given tolerance) to widely used methods for solving general linear systems: the power method and the generalized minimal residual method (GMRES). The numerical results show that one of the proposed algorithms is more competitive that the power and the GMRES methods and over-perform by far the power when \(n\) becomes large.