Tristable and multiple bistable activity in complex random binary networks of two-state units

We study complex networks of stochastic two-state units. Our aim is to model discrete stochastic excitable dynamics with a rest and an excited state. Both states are assumed to possess different waiting time distributions. The rest state is treated as an activation process with an exponentially distributed life time, whereas the latter in the excited state shall have a constant mean which may originate from any distribution. The activation rate of any single unit is determined by its neighbors according to a random complex network structure. In order to treat this problem in an analytical way, we use a heterogeneous mean-field approximation yielding a set of equations general valid for uncorrelated random networks. Based on this derivation we focus on random binary networks where the network is solely comprised of nodes with either of two degrees. The ratio between the two degrees is shown to be a crucial parameter. Dependent on the composition of the network the steady states show the usual transition from disorder to homogeneous ordered bistability as well as new scenarios that include inhomogeneous ordered and disordered bistability as well as tristability. The various steady states differ in their spiking activity expressed by a state dependent spiking rate. Numerical simulations agree with analytic results of the heterogeneous mean-field approximation.


Introduction
Discrete-state stochastic models can be used to describe discrete processes such as the orientation of a spin or the blinking of quantum dots [1]. Additionally, systems with continuously changing dynamics can be mapped to discrete-state descriptions via coarse-graining [2,3,4]. Despite the simplicity of the single units, the collective a Present address: Theory & Bio-Systems Department, Max-Planck-Institute of Colloids and Interfaces, D-14424 Potsdam-Golm, Germany effects of ensembles of coupled units can be highly non-trivial.
In earlier works discrete stochastic three-state models have been used to investigate fluctuation driven spin nucleation on complex networks [5]. Two-and three-state models have been applied to neuronal systems [6,7] and recently to language dynamics [8]. Synchronization behavior, phase transitions and reaction to time delayed feedback [9,10,11,12,13,14,3,4] as well as excitability [15,16] are general aspects that apply to a wide range of natural phenomena.
Most of the referenced works consider Markovian -thereby memoryless -discrete-state models [9,10,11,17]. A disordered environment [18] or the reduction of models with a high number of discrete states to a model with fewer states generically demands a non-Markovian description. Therefore, in continuation of previous work [12,7,15] in this paper a semi-Markovian model [19] of stochastic two-state units is considered. As a new point of interest these units are embedded in an uncorrelated random network whose nodes possess different but independent degrees, in particular this excludes networks with high asortativity or dissortativity. The structure of the network is given by the node distribution p(k). A big number of nodes is assumed, it is known that finite size effects [17] have a strong effect on the dynamic of the stochastic process as well as on the network influence. Complex networks are of top interest in statistical physics because it allows deviation from global coupling without specification of a spatial structure as well as providing a framework to map complex spatial structures to an abstract space. Furthermore, network structures are present in many situations of everyday life, for example transport [20,21,22] and supply networks [23] to mention just a few.
The paper is structured as follows. In Section 2 the master equations for stochastic twostate units are derived. In Section 2.3 the heterogeneous mean-field approximation is used to reduce the set of equations. Section 3 applies the formalism to random binary networks and consists of three subsections. In the first subsection the implicit equations for the steady states are derived and analyzed for saddle-node bifurcations. From this, the scaling of the critical coupling strength due to the network embedding is revealed. The second subsection addresses the limit of vanishing noise. Finally, the last subsection shows the results of the numerical solutions of the mean-field equations in presence of finite noise. Apart from the expected homogeneous ordered bistability that is known from globally coupled units [12,7], inhomogeneous ordered and disordered states as well as tristable states are uncovered. These findings are confirmed via microscopic simulations using the original network structure. In these simulations the network exhibits a variable firing activity by approaching different steady states in dependence on the initial conditions. But the firing of the populations is not synchronized and thus the mean-field stays constant. The firing differs in the spike generating rates for the two populations and whether these are in the rest or excited states. The spike time statistics of individual units are discussed separately. Given that the excited state possess an exponentially distributed waiting time, the spike trains are always nearly Poissonian. For a sharp-peaked waiting time density with no variance, the spike trains become highly coherent.

Two-state processes on complex networks 2.1 Master equation of coupled two-state units
The stochastic two-state units considered in this work can switch between the states 0 and 1. They do so in a stochastic fashion, governed by the waiting time distributions w 0,i (t) and w 1 (t) in the corresponding states, as explained Fig. 1. Focusing on excitable dynamics, state 0 will be called the resting state and state 1 the excited state. In this way, the transition from rest to excited, i.e. from state 0 to state 1, will be assumed as an activation process. The life time in state 0 is exponentially distributed and the state will be left with the transition rate γ.
The backward transition (relaxation from excited to rest) is governed by the waiting time density w 1 (t) to remain in state 1. So far this density is arbitrary except that its mean value shall be τ and does not depend on any other parameters such as the noise intensity, the size and structure or possible dynamical states of the network. In the simulations two specific choices are made. First, the relaxation is treated as a Markovian rate process with a exponential waiting time distribution The two-state units with the waitingtime distributions w 0,i (t) and w 1 (t) are embedded in a binary random network, as indicated by the dashed lines. In this example are five nodes with degree two and four nodes with degree four. The highlighted unit i is one of the latter, hence its activation waiting-time distribution w 0,i (t) is affected by the activity of four neighbors, as indicated by the dotted ellipse and the dotted arrow.
Oppositely, when modeling excitable dynamics, a sharply peaked distribution is more suitable, e.g.
would yield a constant waiting time without any variance in state 1 modeling a fixed delay [12,15,14]. As shown in the appendix A, the bifurcations of the steady states are not affected by the specific choice of w 1 . It depends on the mean time τ spent in the excited state. Only, setups with possible transition to a non-stationary behavior reflect on the choice of w 1 . The two-state units are located on the N nodes of a complex network. Each stochastic element receives the output from the units to which it is linked in the network by edges. This coupling is mathematically realized by introduction of the adjacency matrix A. The activation rate γ i from state 0 to state 1 of the ith node i ∈ {1, . . . , N } is assumed to depend on a signal function f i (t) which contains the adjacency matrix Therein, s j (t) is the output signal of node j depending on its state. In this work, undirected and unweighted networks are considered, hence the adjacency matrix is symmetric with elements A ij = 1, if the units i and j are connected, otherwise A ij = 0. Typically, excitable systems stay in the resting state where no output is produced. Upon sufficient excitation they drastically change their intrinsic dynamics which is emitted as a signal. Here such "spiking" is modeled by a two-valued output function s j (t), which can take the values 1 in the excited state and 0 otherwise. This setting is motivated by neuronal activity or excitable lasers. A symmetric choice would be more appropriate in order to model magnetic spins.
Eventually, in accordance with previous assumptions, the normalized waiting time distribution density of the activation with the time dependent rate from (4) is given by the exponential function Let P i (0, t) denote the occupation probability of unit i to be in state s i (t) = 0 at time t. Analogously P i (1, t) for state s i (t) = 1. Then the balance of probability flows yields the generalized master equationṡ for all i ∈ {1, . . . , N }, where J 0→1,i (t) gives the probability flow from state 0 to 1 of unit i at time t. Since the transition 0 → 1 is a rate process, its probability flow is simply given by The second probability flow is given by all the probability that has flown into state 1 up to time t and stayed there for a time, which is given by the waiting time distribution w 1 (t). Thus it is the convolution of J 0→1,i (t) and w 1 (t ), (8) Using the normalization condition at a given node the temporal evolution of the occupation probabilities P i (1, t) is then given by (cf. (6)) This is a set of N coupled linear integrodifferential equations for the P i (1, t). The complexity is given by the chosen adjacency matrix A i,j which links the almost N equations by the the signal function, see Eq. (4).
By applying a heterogeneous mean-field approximation as described in Section 2.3, the structure of this set of equations is much simpler and the number of differential equations in the set reduces significantly. As consequence of this approximation, the set will become analytically treatable for special cases in the stationary limit.
But before describing this approximation, an appropriate notation for the values describing the dynamical behavior on a complex network will be introduced.

Master equation on complex networks
To include the network topology into the description, it is assumed that the complex network structure can be described as a random network. Central value in this description are the degrees k i of the N nodes. For a given adjacency matrix A i,j with values 0, 1 the degree k i of the ith node is the number of existing links to other nodes of the network. It becomes Let k min and k max be the minimum and maximum degree occurring in the network, respectively, while N k is the number of units with degree k ∈ [k min , k max ]. In a random network the degrees can be treated as random numbers. Their occupation probability p(k) is defined by N k . Then the occupation probability reads Mathematically, the limit make sense if N k ∝ N . Afterwards, sorting the nodes with coinciding degrees k in the network gives the joint probability that any node in the network has this degree and is, respectively, in the dynamical state s = (0, 1) at time t Using Bayes' theorem we split off the occupation probability p(k) as Normalization reads Hence, for the conditioned probabilities the degree is fixed and it becomes P (0, t|k) + P (1, t|k) = 1.
The chosen conditional probabilities P (1, t|k) neglect that various nodes with the same degree k might be linked to k nodes with different degrees. Therefore, without a further approximation, they are so far not suitable for the description of our situation.
To proceed, the signal function (4) which is coupling the ith node to the other nodes will be considered now. Replacing therein the number of the linked node j by the specific degree k j of this node, it becomes Sorting different degrees gives Every node receives its input depending on the specific linked environment. Due to the disorder contained in the adjacency matrix A i,j , the summation is still in general different for distinct nodes. The reduced description for given degrees requires an additional approximation.

Heterogeneous mean-field approximation
As proposed in [24,25,26], the complex network with the adjacency matrix A i,j and with the given degree distribution p(k) will be replaced by a fully connected network with weighted edges. The latter with adjacency matrixÃ i,j shall possess the same distribution of the degrees at the nodes. In detail, it is required that the degreevalues of all nodes in the new network shall coincide with the corresponding ones in the original network. Accordingly, the following shall hold Whereas the sum at the l.h.s. is running over values 1 and 0 of A i,j , at the r.h.s. the sum goes over rational numbers. The assumption is made that the probability to have an edge between the ith and jth node is proportional to the product of the degrees of these nodes, i.e., to k i k j . This is strictly valid only for uncorrelated networks. Further on, taking into account the conservation of degree as required in (19) the following replacing for a fully connected network is defined Figure 2 illustrates the replacement. As a result the single nodes with given degree k couple uniquely to a mean filed. Hence, in contrast to the original network, the fully connected network will allow a mean-field representation with respect to all edges with coinciding degrees k. In [25] the validity of this replacement procedure was discussed and it was successfully applied to complex networks with continuous phase oscillators at the nodes. The approach [24,25,26] is here generalized to discrete stochastic two state units.

Fig. 2:
Representation of a random binary network with k 1 = 12, k 2 = 6 and N = 25 before and after applying the heterogeneous mean-field approximation.
As a result of the replacement, the signal function in (18) is approximated as (21) It follows immediately, that the signal at the node with degree k i coincides for all nodes with the same degree. The node value itself enters only multiplicatively into this expression. Hence, the denominator of the sum is equal for all nodes.
Crossing from the summation over all nodes to sums with the same degrees similar transformations as above are made. The different degrees k i are again divided into classes of units with the same degree k and with occupation number N k . Obviously, kmax k=kmin N k = N has to be satisfied. The denominator can be rewritten In the limit of large number of nodes N → ∞ this expression becomes N k where the symbol · = k · p(k) assigns averaging over the degree distribution. Hence, the signal function of the ith node becomes Therein, the mean-field amplitude r(t) is defined by It is assumed that, after forgetting initial conditions, units with the same degree share stochastic pulse sequences which are statistically identical. Therefore, the same pulse sequence can be assigned to units of the same degree class by taking the average over the corresponding class Therein, the sum runs over N k items due to the action of the δ-function. Since the number of nodes with degree k scales as N k ∝ N , application of the limit of large N yields which was introduced in (14). Thus, the mean-field r(t) introduced in (24) becomes in this limit (27) Eventually, by inserting (23) and (24) via (27) into the mean-field description (15) and (10) the following expression is obtained This is now a set of coupled non-linear integro-differential equations, since the mean-field r(t) depends on P (1.t|k) via (27). Though the structure of the equations looks similar to the previous Master equation (10) it is qualitatively different. As the result of the replacement of the adjacency matrix, the number of equations has reduced drastically compared to the system (15) and (10). The index k in (28) is only running over the different possible degrees in the network k ∈ [k min , k max ] whereas in (15) and (10) it runs over all nodes N . In addition, the dependence of the activation rate γ on the degree k appears uniquely in the argument for all nodes as linear factor in the signal function, cf. (23). In the master equation for nodes with the degree k it reads

Stationary behavior of excitable units
Depending on the specific structure of γ this equation can be highly non-linear. Here in this manuscript, the activation rate γ is assumed to follow Arrhenius' law [27,28]. The two-state system shall mimic the behavior of stochastic excitable dynamics [2] with state 0 being the rest state and state 1 the excited one, respectively. Transitions to the excited state 1 is achieved by overcoming a threshold with barrier ∆U under the influence of noise with intensity D. The corresponding Arrhenius' law of the rate reads with a constant γ 0 defining the time scale. The coupling between units is assumed to be purely excitatory and thus each coupled unit that is already in the excited state will lower the potential barrier by an amount proportional to the coupling strength σ, which is the same for every unit throughout this paper. The following ansatz for the potential ∆U combines the above information Setting p(k) = δ k,N restores the globally connected network with r(t) → P (1, t|N ), which has been earlier studied in detail [14]. Alternatively, it is possible to consider a discrete number of K degrees k i , i = 1, . . . , K. The corresponding degree distribution follows as Insertion of the specific rate and degree distribution (28) yields a set of K nonlinear master equations. A qualitative discussion of the possible stationary solutions P * (1|k i ) which will be approached as t → ∞ will be outlined. Since the equations are nonlinear there might exist a different number of stationary solutions with different stability. If distributed by (32) and with the rate function (31) these stationary states are defined by the set of K coupled nonlinear algebraic equations with the stationary spiking rates of the elements with degree k i and the stationary mean-field amplitude with (33) inserted and averaged over the discrete distribution (32) The maximal number of possible stable solutions of this set of equations can be estimated to be of the order O(2 K ). The behavior is similar to a spin chain with K elements. Every population is stable in the rest state s i = 0. If being coupled, every population with given degree reaches a stationary probability to be in the excited state.
The precise number depends on the specific degree values, the noise intensity D and the coupling constant σ. Generally for low coupling, respectively for high noise only a single solution exists, which is the disordered state of all populations. Lowering of noise, respectively increasing coupling, enlarges the number with multiple stable states, including inhomogeneous cases and the two homogeneous situations where all populations are in the rest or in the excited states. In general, it depends on the initial conditions which state will be populated.
Corresponding to the selected steady-statesolution the spiking rates differ. The spiking rate from rest to excited is determined by the stationary states of the population which the element belongs to. It is expressed by the rate given in (34) where the specific solution has to be inserted. Such state dependent dynamical behavior of neuronal activity was recently discussed for phase oscillators in [29].
In this simplified model no Hopf-bifurcation can take place. All interactions have an aligning effect of the elements similar to a spin chain. Therefore, the existence of stable oscillating, chimaera state, cluster synchronization or chaotic solutions can be excluded. But adding delayed feedback of the mean-field, mixtures of excitatory and inhibitory acting units of the network or systematic shifts in the signal function might be a source for more complex situations.
Special initial conditions (for example all units in the excited state) cause damped oscillations of the mean-field. Even in the case of a δ-function as waiting time density, the assumed exponentially distributed activation times scatter the individual spins with large dispersion. Hence the coherence of special initial conditions is destroyed after a few excitations yielding a stationary meanfield.
Nevertheless, as will be seen in the next section, the individual spins can fire with a small CV 1 resembling oscillatory behavior. In the states with high mean-field values the activation time becomes negligible small. In these situations the spiking is dominated by the recovery time from the excited to the rest state. If this time does not possess remarkable dispersion the CV become vanishingly smalls.
In the next Section 3 further insight into the consequences of several degrees will be given by dealing with the simplest case of network with two populations with different degrees k 1 , k 2 , only.

Random binary networks
In the following, the properties of a binary random network [30,31,26] will be studied. It allows a full sketch of the possible bifurcation scenario appearing in these networks. These ones are randomly connected with two different degrees k 1 and k 2 .
In terms of the degree distribution for a binary network it is supposed that where ν ∈]0, 1[ is the fraction of nodes with degree k 1 . In this paper k 1 > k 2 is set, but due to the symmetry (k 1 , ν) ←→ (k 2 , 1 − ν) the alternative case is also included. The set of two master equations for binary networks read: where we have denoted γ k1 (t) = γ[k 1 /N r(t)] and γ k2 (t), respectively. Equations (37) can be brought into the integral form [32]: supplemented by initial conditions. Therein, z 1 (t) is the survival probability of state 1, Equations (38) have to be supplemented by initial conditions.

Qualitative discussion of steady states
Equations (38) are suitable for calculating the steady states of this system. For a steady state lim t→∞ P (1, t|k i ) = P * (1|k i ) applies. Using (38) and integration by parts gives the following coupled implicit equations for the steady states The values of P * (1|k 1 ) * and P * (1|k 2 ) * define the stationary order in the two subpopulations. In Eq. (40) we introduced the mean relaxation time of the excited state τ = ∞ 0 t w 1 (t)dt and the steady state activation rate γ * ki = γ[k i r * ] for i = 1, 2 depending on the steady state mean-field r * . Equation (27) defines the order parameter of the full network. It becomes r * = kP * (1|k) / k . Taken at steady state, this yields a transcendent equation for the steady state value r * of the mean-field which yields This equation can possess several solutions which we will discuss in detail, later on. In case of large noise D → ∞, only the homogeneous disordered solution r * = 1/2 exists. In this limit the exponential function becomes unity and since we will select γ 0 τ ≈ 1 the disordered state is characterized by r * = 1/2. A constant value of r * does not imply that the activity of the individual nodes has ceased. It is rather the mean activity or flow that is constant (cf. Fig. 8). Oscillating behavior of r * would correspond to synchronization among the units. But without further ingredients like delayed feedback or additional inhibitory coupled nodes such states cannot be reached. The transition to the disordered state can be studied in more detail. Demanding that the first derivatives of l.h.s. and r.h.s. with respect to r coincide at r = r * , provides a condition for a saddle-node bifurcation. The homogeneous disordered state becomes unstable and two new stable solutions occur.
Execution of the derivatives in (41) results in Using (40) this becomes Changing the variables to x k1 := P * (1|k 1 ) − 1 2 and x k2 := P * (1|k 2 ) − 1 2 and rearranging the equation gives Equation (45) defines an ellipse with semi-axes and similarly a 2 with the substitutions k 1 → k 2 and ν → (1 − ν). The ellipse reduces to a point where the two bifurcations merge. It corresponds to P * (1|k 1 ) = P * (1|k 2 ) = r * = 1 2 , at For the γ k given by (31), (47) results in Comparing this to the well known result of all-to-all coupled stochastic two-state units [14], it is visible that it differs only by a scaling factor given by the ratio of the first two moments of the degree (density) distribution. This is a typical network effect in mean-field coupled oscillators [33,34,35,36]. The factor can be interpreted as the mean of the degree distribution of the nearest neighbors [37] assuming that the local structure of the network is treelike. Therefore the effective coupling strength will be introduced. Note that it is evident from (47) that this scaling is only obtained if γ depends exponentially on the mean-field r. In case of low noise a more detailed picture with possibly multiple solutions and ordered states occur. These solutions can be discussed solving (41) graphically and plotting the r.h.s. versus the l.h.s as presented in Fig. 3 for typical situations.
The l.h.s. of equation (41) is a straight line and unbounded whereas the r.h.s. grows monotonically and is bounded between values of the interval [0, 1]. Hence solutions r * are also in this interval. Solutions with one, three or five intersections can be found. Bifurcations between these monostable, bistable of tristable behavior are saddlenode bifurcations or, if these coincide, a pitchfork bifurcation.
If the number of degrees would be increased, the number of steps will also increase in the same manner, giving rise to even higher multistable states.

Solutions with vanishing noise
It is illustrative to look first in detail at the case of vanishing noise. Then the r.h.s. of equation (41) vanishes ∝ exp(−1/D) as r → 0 and approaches unity for large values of r. In between, the r.h.s. makes two jumps with magnitude νk 1 / k and (1 − ν)k 2 / k , respectively. These steps are located at r 1 = 1/(σk 1 ) and r 2 = 1/(σk 2 ). For the r.h.s. to possess one intersection (monostability) with the straight line r, we obtain the following conditions by using that k 1 > k 2 : 1 For vanishing noise the monostable state is always the ordered one with r * = 0, i.e. neither of the two populations is excited.. If one of these inequalities is violated, the mean-field dynamics exhibit bistability. Given that the first one does not hold, besides the ordered solution with r * = 0, a second stable inhomogeneous state appears. The higher degree population is in the excited state and the lower degree population remains in the non-excited one (cf. Fig. 4b). In contrast, if the second inequality is violated bistability occurs between the homogeneous non-excited states and the homogeneous excited ones (cf. Fig. 4a). Finally, if both inequalities do not hold, the solution has five intersections according to a tristable solution between the two homogeneous situations and the one non-homogeneous one (cf. Fig. 4c).

Solutions with finite noise and simulations
Examples of the qualitative behavior of the steady states for various noise levels D are presented in Figs. 4a-4c. The graphs have been obtained by numerical solution of (41) and (40). The main difference of these graphs is the ratio α = k 1 /k 2 . In Fig. 4a, α equals 3/2 and for high noise the disordered state with mean-field r * = 1/2 is stable. The latter bifurcates for lower noise to the known bistability of ordered states [14] which approach (0, 0), (1, 1) as D → 0. These have been written in terms of the state vector P 1 (t) = (P (1, t|k 1 ), P (1, t|k 2 )) T . With respect to the network these solutions are homogeneous states since both populations of the network are ordered in the same states. Figure 4b presents the qualitative behavior with a strong mismatch of the degrees of the two populations, namely α = 4. The graph shows that in this parameter region bistability of the two ordered states occurs for lower noise values. With vanishing D-values the states become (0, 0) and (0, 1). In difference to the previous case, the second solution (0, 1) is inhomogeneous with respect to the two populations in the network. In the (0, 1) state one population is ordered in state 0 whereas the other approaches an ordered state with mean activity 1. It is a result of the strong mismatch α of the degrees and of the asymmetry ν of the two populations. If, for example, the first smaller population with a higher degree is ordered in the excited state 1, it is not able to excite the second larger population anymore. The latter remains in the ordered rest state 0.
Also the coexistence of both scenarios is possible and give rise to a tristable parameter region as presented in Fig. 4c. Here a moderate value of α = 2 was selected. Lowering the noise intensity, three different regions are visible. First, for high noise D > 0.5 the monostable disordered solution exists which is apparent in all figures. The typical bifurcation into two homogeneous states. This is also observed in networks with all-to-all topology. b) k 2 /N = 1 8 : This bifurcation is similar to Fig. 4a but the lower degree population cannot be in the excited state anymore. c) k 2 /N = 1 4 : In this bifurcation diagram there is a monostable regime for large D, a bistable regime with a total ordered and a partial ordered state for intermediate D and the tristable state with three different ordered states for low D.
Lowering the noise this state becomes bistable between the ordered homogeneous states where both populations are in state 0 and an inhomogeneous network where the population that is smaller and stronger connected is with high probability in the ordered state 1 but the larger less connected population is still disordered. This is another type of bistability, this time between an homogeneous ordered state and an inhomogeneous disordered state. Finally, by decreasing the noise intensity D further, the third region is entered and the solutions become tristable between (0, 0), (1, 0) and (1,1) in case of vanishing noise. Stability of the state (1, 0.5) as visible in the intermediate region of Fig. 4c is an interesting event, because it means that the population of the network with the higher degree is in an ordered state whereas the population with the lower degree is in a disordered state. Such partly ordered states have also been reported for the Ising model on correlated scale-free networks [38]. These should not be confused with a chimera state, because the units of the subpopulations are not identical and the value of r does not reflect the synchrony of the phases among the units.
In Fig. 5 the distribution of the different stability regimes for varying degree mismatch α and relative concentration ν are shown for D = 0.1. The tristable region forms an island surrounded by the different types of bistability and connected to the monostable shore by a very narrow region. Going around the island in an anti-clockwise manner one starts at the homogeneous ordered configuration which gradually becomes more disordered and inhomogeneous with maximal disorder in the middle of the right hand side. After the turning point it gets ordered again but this time in the inhomogeneous regime.
The presented findings have been confirmed by microscopic simulations of the coupled network, see Eqs. (6) with w 0 (t) from (5). Numerical investigation of the random binary network is done by solving the Master equations (37) in the Markovian case, namely w 1 (t) as in (1). As shown in (42) and (40), the equations depend on the first moment of the waiting time distribution rather than the shape of the distribution, although higher moments may play a role for other bifurcations (cf. Appendix A). In addition microscopic simulations of a random binary network with 6000 nodes and full network topology, which corresponds to (10), confirm the made approximations.
Exemplary results are shown in Fig. 6 -7 which reproduces Fig. 4c with two different waiting time distributions. The exponential waiting time distribution which was also used in Fig. 4c and a δ-distribution with same mean but no variance. Fig. 8 shows the activity of arbitrary nodes in the two populations in the three different regimes of stable states. The spiking activity of the nodes is presented as symbols versus running time accordingly to the noted degree values k 1 , k 2 at the r.h.s of the graph.
The firing activity in the two different states is different. As discussed already earlier in Section 2.4, units fire seldom in the rest state, but rapidly in the excited state. This dynamical behavior survives also in the inhomogeneous case. In Fig. 8a we present the activity of the units, with an exponential relaxation waiting time distribution. High disorder of the spiking activity in the excited states is the consequence. The CV of the simulated activity is close to 1, which is the value for an independent Poissonian spike train. Differently, in Fig. 8b spiking events from simulations with a δ-distribution are shown. For states where a population i = 1, 2 is in the excited states, i.e. if P (1|k i ) ≈ 1, the measured CV possesses values close to zero . This corresponds to a perfectly oscillatory behavior of the units. The period of this spiking coincides with the time τ the unit stays in the excited state. After this period the unit flips to the rest state. The exponentially distributed time to flip back in the excited state vanishes and also its variance. In consequence, the units behave loike oscillators.
It is important to stress that the constant mean-field value does not correspond to synchronization of the individual units in the excited states, even in case when they practically oscillate. In the rest state the measured CV for both choices of waiting time densities are close to 1.

Conclusion
In this paper, we have investigated semi-Markovian stochastic two-state units embedded in a complex network. A theoretical framework has been developed through a heterogeneous mean-field approximation, which is valid for random uncorrelated networks. Our work thus represents an extension of previous studies on globally coupled two-state systems (especially [14], but also [17,13,6]) to two-state systems that have a complex coupling structure.
As an example, we have focused on a random binary network. Thereby we have discovered qualitative changes in the behavior of the steady states. Specifically, structurally new conformations have been found, including tristable and partially ordered states. Additionally the influence of the network on the critical coupling strength has been revealed. We have corroborated all our theoretical results via numerical simulations. We found that the "Markovianity" of the underlying process has no great effect on the positions and the number of steady states and their bifurcations in this simple setting, but still their basins of attraction may be different. Instead the first moment of the waiting time distribution has the greatest impact and higher moments occur only in bifurcations that have at least co-dimension two. It remains for future studies to pursue our analysis explicitly in cases, where more than two different degrees exist. It will be particularly interesting to see how our findings regarding the multistability will generalize.
Networks of stochastic two-state units can be seen as a toy model for magnetic spins, neurons, blinking phenomena or two valued opinions. The occurrence of tristability is also reported in molecular switches [39,40] and in systems of polaritons [41,42] giving hope to perform ternary logical operations in the future. Hence, we expect that our results are relevant for these real-world systems, where the model considered here may serve as an idealized version. The analytic tractability is a strength of our system, but we believe that still many important extensions await consideration, e.g. including network correlations, more sophisticated waiting time distributions or coupling functions.
As underlined previously, the main assumption behind the heterogeneous mean-field approximation is the lack of degree correlations. Random binary networks tend to be disassortative, i.e. nodes with different degrees are preferentially connected as discussed in [31]. However, for the network examples considered here and for the chosen parameters, these degree correlations become negligible. Therefore, the key assumption is not violated in our study which gave the reason for our analysis. Extending our theory towards correlated networks remains a challenging open problem.   0). a) Network with exponential w 1 corresponding to Fig. 6. b) Network with δ-distributed w 1 corresponding to Fig. 7.

A Derivation of the characteristic equation
To study the stability of steady states by a characteristic equation, we introduce the vector P 1 (λ) which is the Laplace transform of the time dependent deviations δP 1 (t) = P 1 (t) − P * 1 at a steady state. Then, the linearized version of (37) reads in Laplace space λP 1 (λ) − P 1 (t = 0) = J(λ)P 1 (λ), with the Jacobian J(λ). Its formal solution is The final value theorem can be used to calculate the steady states of this system with det(M) = λ 2 − λ(1 − w 1 (λ)) The final value theorem states that the limit lim λ→0 λP 1 (λ) is unique if and only if the denominator of P 1 (λ) has roots with negative real parts and not more than one pole at the origin. Thus indicating bifurcations when one (or several) roots of det(M) cross the imaginary axis. Therefore det(M) = 0 is called the characteristic equation. The fact that w 1 (λ) is the moment generating function of w 1 (t) means that w 1 (λ) = ∞ k=0 t k k! λ k where t k is the k-th moment of w 1 (t). Two typical examples are the pairs Since t 0 = 1, the term (1 − w 1 (λ)) = ∞ k=1 t k k! λ k and thus adj(M) is of first order in λ. Given this information it is clear that P 1 (λ) has only one pole of order one at the origin.

B Alternative derivation of (42) using the characteristic equation
To look for saddle-node bifurcations the lowest terms in λ of equation (55) will be collected. Identifying t = τ results in as the condition for a saddle-node bifurcation. Application of the chain rule in the derivatives, i.e. ∂γ/∂P = ∂γ/∂r ∂r/∂P finally yields This is indeed the same equation as (42). This derivation has the positive side effect that with the aid of the characteristic equation all other bifurcation scenarios can be investigated as well.