The phase diagram of random threshold networks
Threshold networks are used as models for neural or gene regulatory networks. They show a rich dynamical behaviour with a transition between a frozen and a chaotic phase. We investigate the phase diagram of randomly connected threshold networks with real-valued thresholds and a fixed number of inputs per node. The nodes are updated according to the same rules as in a model of the cell-cycle network of Saccharomyces cereviseae [PNAS 101, 4781 (2004)]. Using the annealed approximation, we derive expressions for the time evolution of the proportion of nodes in the “on” and “off” state, and for the sensitivity . The results are compared with simulations of quenched networks. We find that for integer values of the simulations show marked deviations from the annealed approximation even for large networks. This can be attributed to the particular choice of the updating rule.
Threshold networks can be used to model gene regulatory networks Bornholdt and Sneppen (2000); Li et al. (2004); Sevim and Rikvold (2008). The nodes of the network represent genes, and the directed links between them represent interactions between genes. Each node can be in two different states (“on”, “off”). That means that the gene is either expressed or not expressed. Furthermore, each node receives inputs from randomly chosen other nodes that regulate its activity cooperatively. The interactions between the nodes can be excitatory or inhibitory so that one node can activate or repress the expression of another node. In this paper, we study a threshold network where the time development of the states of the nodes is given by the following equation
Here, is a threshold that is the same for every node. The couplings are with equal probability, if node receives no input from node . The input from node to node can therefore take three different values: or . A node becomes activated when the sum of its inputs exceeds the threshold value, and it becomes inactive when the sum of its inputs is below the threshold. When the sum of the inputs gives exactly the threshold value the node does not change its state in the next time step. The nodes are updated in parallel. These dynamics with and a value that varies from node to node were used to model the cell-cycle network of Saccharomyces cereviseae. This model was able to reproduce the overall dynamic properties of the real network Li et al. (2004). There exist several variants of threshold models. In other variants, the can be continuous quantities chosen at random from some probability distribution; the spin values may be instead of 1 and 0 (see for instance Kürten (1988); Rohlf and Bornholdt (2002); Rohlf (2007)), or the update rule in the case that the sum of the inputs is exactly at the threshold can be different.
Models that use spin values can be mapped onto models with spin values by making the substitution . For our update rule (1), this leads to
This means that each node obtains its own threshold value , which depends on the values of the . Therefore the dynamics of the model studied in this paper is different from that of the model studied more widely.
Similarly to random Boolean networks, random threshold networks show a transition between a frozen and a chaotic phase when the network parameters are varied. In the frozen phase, a perturbation at one node propagates during one time step on an average to less than one other node. In the chaotic phase, the difference between two initially almost identical states increases exponentially fast, because a perturbation propagates on an average to more than one node during one time step. In the frozen phase, the length of attractors (i.e., the number of states on attractors) is either 1 or very small. Most of the nodes are frozen, that is they do not change their states anymore in the stationary state. In the chaotic phase, attractors are very long on average, and a non-vanishing proportion of the nodes change their states on the attractors. This phase transition was previously studied in threshold networks in Kürten (1988); Rohlf and Bornholdt (2002); Rohlf (2007).
Ii The phase diagram
With the help of the annealed approximation introduced by Derrida and Pomeau Derrida and Pomeau (1986), one can determine the parameter values and for which the networks are in the chaotic or in the frozen regime. This approximation neglects that the input connections to nodes are constant in time (quenched). It describes therefore a situation where the connections are changed randomly in each time step. The annealed approximation also neglects fluctuations and can therefore become exact only for infinitely large networks (if at all). The parameter , called sensitivity Luque and Solé (1997); Shmulevich and Kauffman (2004), which is times the average probability that the output of a node changes when one of its inputs changes, discriminates between the two phases. If , the network ensemble is said to be in the frozen phase. If it is in the chaotic phase. For the networks are critical. In order to determine , one has to know , the proportion of nodes in state 1 at the considered moment in time. is a function of and becomes constant only when has reached a fixed point.
The annealed approximation has been used successfully to predict the phase diagram of various classes of random Boolean networks. In those networks, correlations between nodes are apparently irrelevant for the evaluation of and . We will see further below that this is not correct for all threshold networks.
ii.1 Time evolution of
Let us first calculate as function of using the annealed approximation. For non-integer , the value of a node will be 1 in the next time step if the sum of its inputs is larger than . Therefore,
Here, is the number of input nodes with value 1. For positive threshold values at least of the input nodes must be active if the sum of the inputs shall be larger than the threshold. For the same reason, the number of positive (excitatory) couplings from these active input nodes must be at least . For negative values all configurations with also contribute to the sum. There are different possibilities to choose active nodes among the input nodes and different possibilities to choose excitatory links among the links from these active nodes. Finally is the probability that input nodes are in state and the others in state , and is the probability that positive and negative couplings are distributed as they are.
For integer-valued , the sum of the inputs can be exactly at the threshold, which is not possible for non-integer . Within the annealed approximation, a node with the inputs at the threshold will be “on” with a probability in the next time step. Equation (3) obtains therefore a second term when is integer and becomes
Here, is again the number of active input nodes with positive couplings. The number of active nodes with negative couplings has to be in order to place the sum of the inputs at the threshold.
Having established the recursion relation for , one can plot maps vs. for different and . The fixed points of the map (see figure 1) are stationary solutions of the annealed approximation.
A fixed point exists whenever the smallest contributing to the sum in (3) is larger than 0. This is the case for all . The fixed point is stable when the slope of the map at is smaller than 1, which is the case for all . For , the map is to leading order in , and the fixed point is therefore unstable for . For , we have to include the next order in , which gives , and therefore the fixed point is stable. For , we have to leading order , and the fixed point is therefore unstable.
In order to obtain information about other fixed points, we iterated numerically the recursion relations for and plotted the map. We found that for and sufficiently small , the only stable fixed point is , but for growing a second stable fixed point appears. Figure 2 shows how the map changes with increasing when . A second stable fixed point with appears at . It moves with increasing slowly towards the value , which is the asymptotic value for . For and , the only fixed point is . For , this fixed point is unstable, as mentioned before, and there exists a stable fixed point with a value . It moves towards with increasing . For , there is only one stable fixed point for all values of . For , the stable fixed point lies between 0.5 and 1 and moves towards 0.5 with increasing .
For , all fixed points can be determined analytically. For , the input of no node can be above the threshold, and there is no fixed point besides . Similarly, for the inputs of all nodes are above the threshold, and therefore . Evaluation of the recursion relation for the remaining values of gives for and for and for and for .
Having found a fixed point for a pair of parameter values and , one can determine the corresponding value of . For non-integer , the change of an input can affect the output only when the other inputs sum up to be directly above or directly underneath the threshold. This leads to
Here, is again the number of active input nodes. The number of active input nodes with positive couplings, , is chosen in such a way that is close to the threshold (directly underneath or directly above). The factor is the probability that the th coupling has the proper sign.
For integer , the situation is again different. A change in an input can affect the output only if the sum of all inputs was before the change. In the opposite situation, where a change in an input places the total input exactly at the threshold, the output does not change. (If we took the annealed approximation to its extremes and ignored the fact that there is a correlation between the state of a node and the state of its inputs, we would need to consider also the case that the sum of the inputs is before the change of one input node.)
We therefore obtain for integer
Using the last two equations, one can evaluate for every combination of and . The resulting phase diagram is shown in Figure 3. Only networks with and a threshold value are critical. Networks with are frozen in the range shown. Where there are two stable fixed points for networks are frozen at and chaotic at . For , networks with integer are more ordered than those with non-integer .
Iii Numerical Simulations
We performed computer simulations of quenched threshold networks for different and and compared the results to those obtained in the framework of the annealed approximation. The number of nodes was in all simulations. The threshold values were chosen in the range , with non-integer thresholds chosen to be . All non-integer thresholds that lie between the same two integers lead obviously to the same dynamical behaviour, therefore it is sufficient to consider these values.
We will first look at the value found after a sufficiently long time when starting with some initial proportion of 1s in the system.
iii.1 The proportion of 0s and 1s
For non-integer , the simulations are in good agreement with the predictions of the annealed approximation.
: The only fixed point for , , is weakly stable because the map has a slope of 1 at this point. Most simulated networks do not reach this fixed point but run into attractors of varying length with a small of the order of . This is due to the fact that the iteration formula can be applied only as long as is larger than of the order . The negative quadratic term in this equation describes the repressive effect of a second active input with a negative coupling. The decrease of comes to a halt when has become so small that there are no more nodes with two active inputs, which will happen for ever smaller values of when the system size is made larger. The discrepancy between the simulations and the annealed approximation is thus clearly a finite-size effect.
For and , the mean values of the quenched networks are in good agreement with the calculated values (we checked for agreement in three decimal places).
: For non-integer , the nonzero fixed point value appears at values that are in accordance with the annealed approximation. For , this happens at . The average value obtained from our simulations shows however a small deviation of about 1% from the calculated value. We can ascribe this small discrepancy again to finite-size effects, since the slope of the map near the fixed point is close to 1 for the value where this fixed point occurs first. For and , the values are again in good agreement.
: For , a stable fixed point value appears at , just as predicted by the calculations. The value of obtained from the simulations and averaged over the attractor agrees with the one obtained from the annealed approximation in three decimal places. Such a good agreement is also found for and 42.
The case of integer is special, and we will see that in this situation the simulations are not in good agreement with the annealed approximation.
: For networks with a value ranging from 2 to 11, the annealed approximation predicts a single stable fixed point . In contrast, our simulations show that already for there exist stationary states with a larger number of active nodes, and the value approached for large times depends on the initial value (see figure 4). Each curve in figure 4 corresponds to one network realization. Each point in the figure is a fixed point of the dynamics, that is an attractor of length one. As one can see, two different networks with the same values of and show approximately the same behaviour. This means that the function does not depend on the detailed realization of the networks. For networks with , the nonzero fixed point appears when the initial proportion of 1s is around 50%.
Only for , the annealed approximation predicts a stable fixed point . For such values, our simulations give a value that is independent of if is not too small. If one compares the value obtained by iterating formula 4 with the value obtained by averaging over several simulated networks, one finds that the quenched networks have around 41% more active nodes than predicted by the annealed approximation. in the publications cited above, where all nodes have the same threshold value. With increasing , this deviation from the annealed approximation decreases. It is about 25% for and about 20% for .
In order to understand how a broad set of dynamical fixed points can emerge in these systems, we note the following: (a) At every fixed point, there are nodes the sum of whose inputs is at the threshold. If changing the state of such a node does not change the state of any node influenced by it, we have found another fixed point with a different number of active nodes. (b) All fixed points of a network with or with are also fixed points of this network if . One can expect that all observed fixed points should have values between these two boundary values, and this is indeed observed. (c) If there exists a set of fixed points with different values, the final value of reached in a simulation should depend on the initial value. This in turn means that the state of a node at a fixed point depends on the dynamical history of this node and its inputs. Such correlations between the states of a node at different moments in time are not taken into account in an annealed approximation, which is therefore no good approximation in this situation. With increasing , the curves in figure 4 become increasingly independent of . This can be attributed to increasing transient times, which weaken the “memory” of the initial state of the nodes. (d) There can exist sets of active nodes that are connected by positive directed links among each other in a way that loops are formed. When the sums of all other inputs to the nodes in such loops are zero, they are at the threshold and once activated will keep their value. Correlations of this type between nodes are not included in annealed models.
and : For networks with a threshold value , a second stable fixed point value appears at in the annealed approximation, but can be found in quenched networks already for . At , the second stable fixed point is about 36% higher than the value predicted by the annealed approximation, and this difference decreases to 26% and 21% for and . For , a second stable fixed point appears also at a lower (at ) than expected from the annealed approximation, where it appears at . The value of reached at is about 29% higher than the expected value. We find a deviation of about 23% and 20% for and . All observed values are between those obtained for and if , and between those obtained for and if .
For and ranging from to 5 we again find many dynamical fixed points that have a different from the value calculated with the help of the annealed approximation, see figure 5. With increasing however, approaches the expected value.
In Figure 5, one can observe another interesting effect, which occurs also for all other integer and for not too large values. Depending on the sum of and , the curves exhibit different behaviour at high values. If is even, increases at the end as in this case it is more likely that the inputs of a node are at the threshold than for odd .
Considering the trivial case is also instructive: Every node has exactly one input. Starting from a randomly chosen node, we can follow the chain of inputs preceding this node. For large system sizes, the average length of such chains is long (it is of the order of ) Flyvbjerg and Kjær (1988). Every chain eventually ends in a loop. Along the chain, positive and negative couplings follow in a random order. We can easily find the fixed points of such a system of chains: Obviously, having all nodes in state 0 is a fixed point of the dynamics. If we then switch on a node that has only negative output links, we obtain another fixed point. If we switch on a node that has a positive output link, the node influenced through this link must also be switched on, and so on, until the end of the sequence of positive couplings is reached. The fixed point with the maximum number of “on” nodes is obtained by assigning a 1 to all nodes with a positive input link and to all those nodes with a negative input link that are preceded by an odd number of nodes with a negative input link. In this way, all nodes that have a negative input link and an input node in state 1, are in state 0. The value associated with this fixed point is
If we had , this would be the only fixed point. The example thus demonstrates that the maximum possible fixed point value for for integer is identical to the one obtained by slightly lowering the value of . Similarly, the minimum value 0 is the fixed point value obtained if is slightly larger than 0. All intermediate values of for are obtained by switching off part of the “on” nodes in the state with the maximum value.
This example demonstrates also that not all fixed point values of can be reached from a random initial state. For instance, the maximum value 2/3 of cannot be reached by starting from such a random initial configuration. If initially all nodes are in state 1, i.e. if , we have after 1 time step , where all nodes with a positive input link are in state 1. If there are less 1s initially, there cannot be more 1s in the final state. Values between and can therefore only be reached by starting from specially prepared initial states. This observation for explains why the simulations for larger (and also for other integer ) give values between those obtained with the annealed approximation for the neighbouring non-integer values, but not the entire interval between these boundaries is reached from a random initial state with fixed .
The fact that the agreement between the annealed approximation and the quenched networks becomes better for larger can be ascribed to the narrowing of the interval between the boundaries with growing .
Figure 6 shows the final value as a function of the initial value for networks with and different negative . (For the parameter values shown all networks have , that is they are frozen.) As one can see, for non-integer the final proportion of 1s in the network does not depend on . Again the values of obtained with the simulations agree well with those obtained using the annealed approximation. For integer , we find dynamical fixed points with a value that depends on . The mean value of is always smaller than the value predicted by the annealed approximation. The values of for the networks with integer always lie between those of networks that have neighbouring non-integer values with the same , as we have also observed for non-negative integer values of .
iii.2 The phase diagram
Next, we report on the dynamical properties of the simulated networks, that is the lengths of their attractors and the number of nodes that change their states while the network is moving through the attractor. We consider the networks at the different possible fixed point values . A fixed point of in the annealed approximation does not necessarily imply that the network dynamics reaches a fixed point in state space (i.e. a dynamical fixed point). In our simulations, all attractors in state space have a constant value of (with some fluctuations around it because the system size is finite), which means that the proportion of nodes changing their state from 1 to 0 is at each step approximately equal to the proportion of nodes changing their state from 0 to 1, as suggested by the annealed approximation. We compare the results of the simulations with the values calculated within the annealed approximation.
The only parameter combination for which the considered networks can be critical is and . The attractors found in the simulations have lengths of one to four digits, the number of nodes that change their state on the attractors is of the order of , which is compatible with the expected number of the order of Kaufman et al. (2005).
For and , the networks should be chaotic according to the annealed approximation. In the simulations, the attractors are longer than our search range (max. transient length: , max. attractor length: ) which is consistent with the expectation of having a chaotic network. The annealed approximation predicts furthermore that networks with larger are frozen when the only stable fixed point is . Networks with values for which a second stable fixed point exists have at this fixed point and should therefore be chaotic. Simulations of networks with non-integer and 2.5 show results that are consistent with this prediction. Quenched networks with integer show again deviations. For and , even networks with values at the second fixed point are frozen and not chaotic. For and , no attractors are found within the search range, pointing at chaotic dynamics. For the situation is different. As stated in the previous section, a second stable fixed point appears already at , and the networks show chaotic behaviour at this fixed point. According to the annealed approximation, the first chaotic network should have . The same is true for networks with : they are chaotic when they reach the second stable fixed point , but this fixed point appears for values smaller than predicted by the annealed approximation.
For , the situation is similar to that for . According to the annealed approximation networks should be chaotic from on. But in all simulations of networks with connectivities up to we find only fixed point attractors, which means that these networks are in the frozen phase.
For , we first chose the values of and such that the networks are expected to be in the frozen phase, and we found fixed point attractors in all simulations. Then we took a closer look at parameter values close to the transition between ordered and chaotic dynamics (cf. figure 3). For integer , we found again deviations from the annealed approximation. Just as for and 1, the frozen phase is extended to higher values. This means that we find networks with fixed point attractors only in regions of the parameter space where they should be chaotic according to the annealed approximation. For non-integer and , chaotic networks should be found for and respectively. In all simulations with these parameters, the attractor lengths exceeded the search range. For and , on the other side of the phase boundary where networks are expected to be frozen, we find mostly short attractors with less than of the nodes changing their state. The major part of these networks is frozen, in agreement with the calculated value . For and , the situation is similar. For , the annealed approximation predicts chaotic dynamics for . But since for , we are very close to the boundary for this parameter value. The simulated networks have attractor lengths ranging from one digit to values exceeding the search range. For , we have , and the attractors are short, with less than of the nodes changing their values.
To summarize, the phase diagram figure 3 obtained by the annealed approximation is valid for the quenched system for all non-integer . For integer , the transition from the frozen to the chaotic phase occurs at a larger value of than predicted by the annealed approximation, and for integer it occurs at a smaller value. Since most of these transitions do not lie in the window of values shown in figure 3, the corresponding figure obtained from our simulations looks hardly different, therefore we do not include it.
The reason why the transitions from the frozen to the chaotic phase do not occur for the values predicted by the annealed approximation when is integer, is the same as the reason why the stationary values do not agree with the fixed points calculated with the annealed approximation. For integer , the annealed approximation does not capture correctly the dynamical properties of the system, because it neglects memory effects. Even when we evaluate within the annealed approximation by using the values obtained from the simulations, the calculated phase transitions do not occur at the same values as those obtained by computer simulations when is integer.
We investigated the phase diagram of threshold networks with real-valued thresholds and an updating rule that does not change the state of a node when the sum of its inputs gives exactly the threshold value. We compared the analytical results obtained by using the annealed approximation with the results obtained from computer simulations. We evaluated the proportion of 1s in the networks and the sensitivity to changes of the state of an input. We found that the annealed approximation is valid in the case of non-integer thresholds, but that it does not agree with the simulations in the case of integer thresholds. We ascribed this discrepancy to memory effects that are not captured by the annealed approximation. In the studies mentioned before Kürten (1988); Rohlf and Bornholdt (2002); Rohlf (2007), this situation did not occur, and the annealed approximation was sufficient to calculate the phase transitions correctly.
Let us now briefly return to the model of the cell-cycle network of yeast Li et al. (2004), which motivated us to study threshold networks with this special kind of updating rule. This model consists of 11 nodes, and it shows seven fixed points. The dominant fixed point corresponds to the G1 phase of the cell cycle, during which the cell grows. Although the trajectory corresponding to the cell cycle is impressively stable, the dominant fixed point is very sensitive to specific perturbations at certain nodes Fretter and Drossel (2008). When the state of one of the 11 nodes is changed, the network returns to the fixed point only in 6 out of the 11 cases. In the other cases, the dynamics is attracted to one of the other fixed points of the network. 6 of the 7 fixed points can be reached from other fixed points by changing only one node, and there is a group of three nodes amongst which all these changes occur. The inputs of all three nodes are at the threshold, and changing the state of one of these nodes does not change the state of any other node. As we have seen in this paper, such a set of fixed points, which differ by the state of one node, is characteristic of a threshold network with an integer-valued threshold. It is due to the update rule that a node keeps its state when the sum of its inputs is exactly at the threshold. This raises the question whether the non-dominant fixed points have a biological meaning, or whether they are just artefacts of the update rules of the model.
Acknowledgements: This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under contract no. Dr200/4-1.
- S. Bornholdt and K. Sneppen, Proc. R. Soc. Lond. B 267, 2281 (2000).
- F. Li, T. Long, Y. Lu, Q. Ouyang, and C. Tang, Proc. Natl. Acad. Sci. USA 101, 4781 (2004).
- V. Sevim and P. A. Rikvold, J. Theor. Biol. 253, 323 (2008).
- K. E. Kürten, Phys. Lett. A 129, 157 (1988).
- T. Rohlf and S. Bornholdt, Physica A 310, 245 (2002).
- T. Rohlf (2007), arXiv:0707.3621.
- B. Derrida and Y. Pomeau, Europhys. Lett. 1, 45 (1986).
- B. Luque and R. V. Solé, Phys. Rev. E 55, 257 (1997).
- I. Shmulevich and S. A. Kauffman, Phys. Rev. Lett. 93, 048701 (2004).
- H. Flyvbjerg and N. J. Kjær, J. Phys. A 21, 1695 (1988).
- V. Kaufman, T. Mihaljev, and B. Drossel, Phys. Rev. E 72, 046124 (2005).
- C. Fretter and B. Drossel, Eur. Phys. J. B 62, 365 (2008).