Difference equation for tracking perturbations in systems of Boolean nested canalyzing functions

This paper studies the spread of perturbations through networks composed of Boolean functions with special canalyzing properties. Canalyzing functions have the property that at least for one value of one of the inputs the output is fixed, irrespective of the values of the other inputs. In this paper the focus is on partially nested canalyzing functions, in which multiple, but not all inputs have this property in a cascading fashion. They naturally describe many relationships in real networks. For example, in a gene regulatory network, the statement “if gene $A$ is expressed, then gene $B$ is not expressed regardless of the states of other genes” implies that $A$ is canalyzing. On the other hand, the additional statement “if gene $A$ is not expressed, and gene $C$ is expressed, then gene $B$ is automatically expressed; otherwise gene $B$'s state is determined by some other type of rule” implies that gene $B$ is expressed by a partially nested canalyzing function with more than two variables, but with two canalyzing variables. In this paper a difference equation model of the probability that a network node's value is affected by an initial perturbation over time is developed, analyzed, and validated numerically. It is shown that the effect of a perturbation decreases towards zero over time if the Boolean functions are canalyzing in sufficiently many variables. The maximum dynamical impact of a perturbation is shown to be comparable to the average impact for a wide range of values of the average sensitivity of the network. Percolation limits are also explored; these are parameter values which generate a transition of the expected perturbation effect to zero as other parameters are varied, so that the initial perturbation does not scale up with the parameters once the percolation limits are reached.

Figure 1

Hamming distance versus model (left column), for PNCF network with $n=64,d=1$ (first row), $d=2$ (second row), and $d=5$ (third row). The Hamming distance is computed over a number of random initial conditions. Both the network and the model are iterated 100 steps, and the actual corresponding values are plotted in the right column. The initial perturbation is on eight randomly selected nodes. Notice that although there is not a perfect match (due also to rounding errors on small numbers), the model follows pretty closely the network results.

Figure 2

Graphs of the function $f\left(y_k\right)$ given by (3) for $d=10$ and six values of parameter $n$ (see legend). Observe that the graphs are practically superimposed. The red ladder illustrates graphically the first few iterations of (3), generating a cobweb plot.

Figure 3

Graphs of the function $f\left(y_k\right)-y_k$ for $d=10$ and six values of parameter $n$ (see legend). The points at which the graphs cross the horizontal axis represent the second fixed points of the map (3). The crosses indicate the leading order approximation, (8), for $n=20,40$, and 60 (left to right).

Figure 4

Top panels: Lyapunov exponents for $f\left(y_k\right)$ as functions of the canalyzing depth $d$. Bottom panels: bifurcation diagrams for the same values of the parameters. The left panels correspond to $n=10$, while the right panels correspond to $n=50$.

Figure 5

Variation of the dynamical impact $y_k$ across networks of sizes at most $n=20$, with $d=0,1,\cdots ,n$, ranging from no canalyzing through full canalyzation, and initial conditions $y_0$ given by (2). The averages are computed over $T=200$ iterations. The horizontal axis represents the average sensitivities over all $(n,d)$ combinations. We collect statistics of $Q(n,d)$ by values of $s(n,d)$ and plot the means with $10\%$ error bars. Note that the maximum impact is much larger than the average impact for smaller sensitivity, but that they become roughly equal as $\langle f^{\prime }\left(0\right)\rangle$ increases. The right graph is a zoom in on the left graph around the critical sensitivity $\langle f^{\prime }\left(0\right)\rangle =1$. The maximum impact is significantly larger than the average impact around the critical sensitivity $\langle f^{\prime }\left(0\right)\rangle =1$, followed by a clear drop leading to comparable values of the average and maximum impact for larger values of the average sensitivity. Thus the average sensitivity has a significant influence on the dynamical impact.

Figure 6

Analog of Fig. 5 with focus on the variation of the dynamical impact $y_k$ across networks of sizes at most $n=400$, with fixed ratios $d/n$ as specified in the plots, ranging from reduced canalyzation through high canalyzation, and initial conditions $y_0$ given by (2). The horizontal axis represents the average sensitivities over all $(n,d)$ combinations for each fixed ratio $d/n$. We collect statistics of $Q(n,d)$ by values of $s(n,d)$ and plot the means without the $10\%$ error bars of Fig. 5. Note that the maximum impact can be significantly larger than the average impact around $\langle f^{\prime }\left(0\right)\rangle =1$, but that they become roughly the same otherwise. Observe the change in the shape of the graphs as $d/n$ increases, leading to mostly small values of $\langle f^{\prime }\left(0\right)\rangle$ for increased level of canalyzation.

Figure 7

Top: ${\langle y_{500}\rangle }_{y_0}$ averaged over several values of the initial perturbation $y_0$ in the specified interval, versus the canalyzing depth $d$ and various values of the network size $n$ within the bounds specified in the plot. For large enough $n$, all graphs transition to zero around the same value of $d\sim 15$, so the damage does not scale up with increased canalyzing depth. The graphs for $n>4096$ are superimposed over $n=4096$. Middle: Analog of the top figure, generated by replacing $d$ with the ratio $d/n$ which is independent of the network size. The (average) expected damage does not scale up with $d/n$. For canalyzing ratios $d/n>0.35$ the expected damage is basically null for all $n$ values, whereas for large enough values of $n$ the transition to zero occurs for $d/n<0.1$ with approximation. Bottom: Analog of the top figure, but with a switch of $n$ and $d/n$. Notice that for $d/n\ge 0.3$ the transition to zero occurs for very small $n$ values, approximately $n<50$. Even for slightly larger values of $d/n$ the damage still converges to zero for $n<200$.