Nonlinear signaling on biological networks: The role of stochasticity and spectral clustering

Signal transduction within biological cells is governed by networks of interacting proteins. Communication between these proteins is mediated by signaling molecules which bind to receptors and induce stochastic transitions between different conformational states. Signaling is typically a cooperative process which requires the occurrence of multiple binding events so that reaction rates have a nonlinear dependence on the amount of signaling molecule. It is this nonlinearity that endows biological signaling networks with robust switchlike properties which are critical to their biological function. In this study we investigate how the properties of these signaling systems depend on the network architecture. Our main result is that these nonlinear networks exhibit bistability where the network activity can switch between states that correspond to a low and high activity level. We show that this bistable regime emerges at a critical coupling strength that is determined by the spectral structure of the network. In particular, the set of nodes that correspond to large components of the leading eigenvector of the adjacency matrix determines the onset of bistability. Above this transition the eigenvectors of the adjacency matrix determine a hierarchy of clusters, defined by its spectral properties, which are activated sequentially with increasing network activity. We argue further that the onset of bistability occurs either continuously or discontinuously depending upon whether the leading eigenvector is localized or delocalized. Finally, we show that at low network coupling stochastic transitions to the active branch are also driven by the set of nodes that contribute more strongly to the leading eigenvector. However, at high coupling, transitions are insensitive to network structure since the network can be activated by stochastic transitions of a few nodes. Thus this work identifies important features of biological signaling networks that may underlie their biological function.

Figure 1

The network activity $p\left(t\right)$ as a function of time for a range of initial conditions. Here, we use an ER network of size $N=500$ that is above the percolation threshold $\left(q=0.03\right)$, and with fixed excitability $\eta =0.005$. (a) Inactive phase with $s=24$. (b) Bistable phase with $s=27$. (c) Active phase with $s=60$. (d) Density plot of the distribution of steady state probability $p_{\infty }$ for the last 50 time units of a simulation up to $T=100$.

Figure 2

The network activity $p\left(t\right)$ for an ER network below the percolation transition $\left(q=0.0015\right)$. Network size is $N=500$ and excitability is fixed at $\eta =0.005$. (a) Inactive phase with $s=5.6$. (b,c) Active phase with $s=15.3$ and $s=50.1$. (d) Density plot of the distribution of steady state probability $p_{\infty }$ for the last 50 time units of a simulation up to $T=100$.

Figure 3

Mean waiting time $\langle t_w\rangle$ as a function of the connectivity $s$. This is the time an inactive system starting with $p\left(t\right)=0$ will take to reach the threshold $p\left(t_w\right)=p_{th}$. The waiting time is averaged over 1000 independent simulation runs. (a) Average waiting time $\langle t_w\rangle$ for an ER network above the percolation threshold and a threshold $p_{th}=0.25$. (b) Average waiting time $\langle t_w\rangle$ for an ER network below the percolation threshold and a threshold $p_{th}=0.20$.

Figure 4

The network activity $p\left(t\right)$ vs time for a range of initial conditions. Here, we simulate a BA network of size $N=500$, with connectivity $m=5$, and fixed excitability $\eta =0.001$. (a) Inactive phase with $s=20$. (b–d) Bistable regime with $s=43.2,44.2$, and $45.7$, showing stochastic transitions from the inactive to active phase. Note the long dwell times in each state. (e) Time evolutions of system in the active regime $\left(s=99.4\right)$, where all initial conditions evolve to a high activity state. (f) Density plot of the distribution of steady state probability $p_{\infty }$ for the last 100 time units of a simulation up to $T=200$. The gap separating the inactive and active branches is denoted as $\mathrm{\Delta }$.

Figure 5

(a) Average waiting time to transition between inactive and active branches. Blue line denotes the waiting time $\langle t_w\rangle$ to transition from the inactive branch to the active branch. Black line denotes the average lifetime $\langle t_a\rangle$ spent on the active branch. The red line shows the averaged steady state probability $\langle p_{\infty }\rangle$ corresponding to the active and inactive branches. (b) The gap at onset $\mathrm{\Delta }$ vs network parameter $m$. (c) Density plot of the steady state probability $p_{\infty }$ of a BA network with $m=2$. Network size is $N=500$ with excitability $\eta =0.001$.

Figure 6

(a) Steady state solutions to the mean-field equations with $\eta =0.005$ and $\gamma =2$. (b) Mean-field approximation of steady state probability $p_{\infty }$ as a function of $s$ of an ER network above (red) and below percolation threshold (black). (c) Mean-field approximation of steady state probability $p_{\infty }$ as a function of $s$ for a BA network with connectivity $m=5$. In this simulation we use $\eta =0.003$. (d) Plot of the gap $\mathrm{\Delta }$ vs $m$.

Figure 7

(a) Normalized ordered components of the leading eigenvector $\left({\varphi }_{1i}\right)$, and activity level at node $i$ at steady state $\left({\overline{\eta }}_i\right)$, for a BA network with $m=2$ and $s=40$. Only the 50 largest components are shown, and all parameters are the same as in Fig. 5. (b) Network structure representing all nodes with ${\varphi }_{1i}>0.2$ in the center and the high activity sites with ${\overline{\eta }}_i>0.5$ highlighted in red.

Figure 8

(a) Normalized ordered components of the leading eigenvector $\left({\varphi }_{1i}\right)$, and activity level at node $i$ at steady state ${\overline{\eta }}_i$, for an ER network below the percolation threshold. Parameters are identical to those used in Fig. 2. (b) Network structure showing eigenvectors corresponding to the three largest eigenvalues (enclosed in dashed lines). High activity sites with ${\overline{\eta }}_i>0.5$ are highlighted in red. Here, the network structure shows only the nodes that have more than two connections. The remaining weakly connected nodes were not active at steady state and are not shown.

Figure 9

Localization of the leading eigenvector in the ER network. (a) Leading cluster size $\mathrm{\Gamma }=n_1/N$ as a function of the linking probability $q$. $\mathrm{\Gamma }$ is averaged over 80 independent network configurations. (b) Ordered components of the leading eigenvector above $\left(q>q_c\right)$ and below $\left(q<q_c\right)$ percolation threshold for a particular realization of an ER network with $N=500$. (c) Leading cluster size $\mathrm{\Gamma }=n_1/N$ as a function of system size $N$ below the percolation threshold $\left(q=1/2N\right)$. (d) Leading cluster size $\mathrm{\Gamma }=n_1/N$ as a function of system size $N$ above the percolation threshold $\left(q=3/N\right)$. For (c,d) points are computed by averaging over 500 independent network configurations.

Figure 10

Structure of the leading eigenvector in BA networks. (a) Leading cluster size $\mathrm{\Gamma }=n_1/N$ as a function of connectivity $m$. (b) Leading cluster size $\mathrm{\Gamma }=n_1/N$ as a function of system size $N$. Results are averaged over 150 network configurations.

Figure 11

Signal spreading. (a) Plot of the number of active nodes $n\left(t\right)$ as a function of time in the case where the network is inactive at $t=0$. Here, we use a BA network with $m=2$, system size $N=500$, and $\eta =0.001$. Connectivity is fixed at $s=39.8$. (b) Network graph with leading cluster placed in the interior. Here the highlighted nodes correspond to the ten largest $P_i$ computed using $M=2000$ simulations. The threshold for activation is fixed at $n_c=20$. (c) Same as above with $s=159.2$. (d) Highlighted nodes correspond to the top ten nodes active at a threshold of $n_c=1$.

Figure 12

Pearson correlation between leading eigenvector and the distribution of active sites $\vec{P}$ at threshold as a function of the network connectivity. Points were computed using 2000 independent simulation runs. Red line indicates the average steady state $\langle p_{\infty }\rangle$ from the data shown in Fig. 5.

Figure 13

Density plot of the distribution of steady state probability $p_{\infty }$ for the exponent $\gamma =1,2$, and 4. The network is identical to that used in Fig. 4.