### Cell differentiation and attractor dynamics

Boolean networks [6] have proved effective in representing GRN structure and dynamics in many systems, including *Drosophila* development [12, 13], angiogenesis [14], eukaryotic cell dynamics [15], and yeast transcription networks [16]. Each gene in a network is represented as a node whose regulation by other genes is modeled using updating rules as logic functions. An expressed gene is assigned the value true and a non-expressed gene the value false.

A Boolean network with *n* genes has 2^{n} possible states, denoted as *Ŝ*. At each step in the simulation, the next state *ŝ*_{t+1}∈ *Ŝ* is determined by applying each gene's logic function (representing the regulatory interactions) to the current value of the genes in *ŝ*_{
t
}. Let this computation be defined as *ŝ*_{t+1}*← D*(*ŝ*_{
t
}) where *D*(*ŝ*_{
t
}) is the deterministic mapping function that finds the next state of the network given the current state. As the network is executed by repeated applications of *D*(*ŝ*_{
t
}), the state will reach a previously visited state, and thus, since the dynamics are deterministic, enter into an attractor which represents a fixed point of the system. Attractors can be single states, called point attractors, or consist of more than one state that the network continuously transitions between, called cyclic attractors. Let *â* = *D*^{∗}(*ŝ*) be the resulting network attractor state reached when starting at *ŝ* and applying the logic functions until the attractor state *â* is reached.

In this work, cell types are considered attractors in the state space of possible gene expression profiles [17] and cell differentiation is modeled as the process of transitioning from one attractor to another [18].

### Network construction

A random Boolean regulatory network is generated by randomly connecting a varying number of nodes, then instantiating each node with a randomly generated logic function. To replicate networks found in natural systems, we created only networks that operate in the critical domain, rather than ordered or chaotic. Critical networks implement maximal information flow [3] and have the lowest attractor basin entropy [19]. Evidence that GRN's tend to be critical is given in [20]. To generate critical networks, the parameters are set according to *s* = 2*qp*_{
N
}(1 − *p*_{
N
}) where *s* is the sensitivity of the network to perturbations in gene values, *p*_{
N
} is the probability of the output of each Boolean function being 1, and *q* is the count of inputs to each Boolean function [21]. When *s* = 1 a single bit change is on average propagated to one other node and the network is in the critical domain. In an ordered network, *s* < 1 and perturbations tend to die out, while in a chaotic network, *s* > 1 and perturbations tend to grow. In this work *s* was fixed at 1 and *p*_{
N
} was adjusted depending upon the value of *q*.

The attractors of each random Boolean regulatory network are determined, then Markov chain analysis is performed to determine the transition probabilities between all possible states. This allows determination of the MFPTs between each pair of attractors [9]. The MFPTs allow the construction of a graph whose nodes are attractors and weighted edges are the MFPT value between different nodes. Figure 8.a shows a sample graph. MFPT graphs for cellular differentiation are expected to have a small MFPT value for forward edges (moving from less to more specialized cell types), large values for reverse edges, and large values in both directions for transitions between attractors at the same level of tree (level is the number of transitions from the root). In [5] a method was introduced that applied successively higher MFPT thresholds to prune edges from this complete MFPT tree as a means to identify separation among subsets of close attractor states as illustrated in Figure 8(b). The effects of changing the threshold from low to high was proposed as a possible mechanism for cellular differentiation with the low threshold representing pluripotency and the process of raising the threshold as type specialization as attractors become more and more isolated. This model proposes that cells differentiate by actively controlling their sensitivity of expression noise and can account for the observation that terminally differentiated cell states tend to be more stable than pluripotent states.

### Network search

We perform a uniform Monte Carlo search over the space of critical random Boolean networks. For each network we find the attractors and compute the MFPT between all possible attractor pairs (extended from code posted at http://code.google.com/p/pbn-matlab-toolbox[9]). Using the MFPT values, for each type of multistable switch added to the network, we draw a density plot where the x-axis is the forward MFPT and the y-axis is the reverse MFPT (we consider the edge with lower MFPT as forward). The acquired density plots are used to determine the distribution of directional, non-directional, separated and non-separated probability transitions between attractors in each network.

### Network types

We investigated nine types of networks. Approximately 5 ∗ 10^{4} networks of each type were explored to find 5, 000 networks of each type with five or more attractors. The different network types come from the use of the 4 multistable switches that are shown in Figure 4. The first switch (Figure 4 a) is a bistable switch, a small local circuit with feedback loops. This is a common switch in biological networks and it controls binary branch points between two mutually exclusive cell lineages[1, 10]. The truth table of the functions in this switch is [1,1,0,1] for binary numbers [00,01,10,11] respectively. The other three switches all encode mutual inhibition between two genes. The first is *MI00* and is based on the network synthesized in [10]. *MI00* includes two incoherent feedback loops. The final two switches extend mutual inhibition with the addition of one *MI+0* or two *MI++* positive (coherent) feedback loops. These two switches were explored in [22], where it was shown that the positive feedback loops can introduce additional shallow attractor basins in continuous ODE network models.

The four switches were used as described above to construct nine different types of networks: no motif, one BS, one MI00, one MI0+, one MI++ and then four more network classes each with two of the same switch. Note that when a motif defined in Figure 4 is embedded, two nodes of the original RBN are selected randomly, their logic functions replaced and inputs and outputs rewired. For illustration, consider how a MI0+ motif is embedded into a RBN. Starting with a RBN (see Figure 4 (a)), two nodes are selected randomly and their truth tables are changed to [0,1,0,0]. Then, the 2^{0} input of the second node is wired to the output of first node and, conversely, the 2^{0} input of first node is wired to the output of the second node. The small or-gate and not-gate are not considered in wiring, because they were previously considered in the truth tables of their respective nodes.

### Mean first passage time

The first-passage time (*FPT*), also called first hitting time, is the time taken by a stochastic system for the first visit of a specific state. Mathematically, *FPT* is defined as *F*_{
k
}(*ŝ*_{
x
}*, ŝ*_{
y
}): the probability that starting in state \widehat{x}, the first time the system visits a state *ŷ* will be at time *k*. In the case of Boolean networks, time is the path length of state transitions. Considering *p*_{
xy
} as the probability of transition between states *x* and *y*, then *F*_{1}(*ŝ*_{
x
}*, ŝ*_{
y
}) = *p*_{
xy
}. As equation 1 shows, for *k* ≥ 2, *F*_{
k
} is calculated by a recursive iteration over all transitive relations: for all *z* states in the network dynamics, *F*_{
k
}(*ŝ*_{
x
}*, ŝ*_{
y
}) is the probability of a one step transition from state *x* to *z* times the *FPT* from state *z* to *y* in *k −* 1 steps.

{F}_{k}\left({\widehat{s}}_{x},{\widehat{s}}_{y}\right)={\displaystyle \sum _{{\widehat{s}}_{z}\in {\left\{0,1\right\}}^{n},z\ne y}}{p}_{xz}{F}_{k-1}\left({\widehat{s}}_{z},{\widehat{s}}_{y}\right)

(1)

Probabilistically, there are two possibilities to reach state *y* from *x*; either *y* is a deterministic target for *x* and no bit flips occur due to the noise, or an aggregate of bit flips drive the transition from *x* to *y*. So the equation for *p*_{
xy
} can be written as follows.

pxy=\left\{\begin{array}{cc}{\left(1-{p}_{e}\right)}^{n}\hfill & y\leftarrow D\left({\widehat{s}}_{x}\right)\hfill \\ {p}_{e}^{{h}_{xy}}{\left(1-{p}_{e}\right)}^{n-{h}_{xy}}\hfill & y\leftarrow \eta \left({\widehat{s}}_{x},{h}_{xy}\right),{\widehat{s}}_{x}\ne {\widehat{s}}_{y}\hfill \end{array}\right.

(2)

where *d*_{
ij
} is equal to 1 if there is a deterministic transition from *x* to *y* in the network dynamics, otherwise it is 0; *p*_{
e
} is the probability of a single bit flip resulting from noise and *h*_{
xy
} is the Hamming distance between two states; *n* is the total number of nodes in the network.

Although the *FPT* is a valuable measure, the average time it takes to reach state *y* from state *x*, termed Mean First Passage Time (*MFPT*), is of greater interest. MFPT in Boolean networks was introduced by Shmulevich et al [23] and is defined as:

MFPT\left({\widehat{s}}_{x},{\widehat{s}}_{y}\right)={\displaystyle \sum _{k}}k{F}_{k}\left({\widehat{s}}_{x},{\widehat{s}}_{y}\right)

(3)

A low *MFPT* between two states indicates that starting from the first state, the second state is easily reached by gene expression noise. Figure 9 shows *F*_{
k
}*, kF*_{
k
}, and MFPT for the transition between two arbitrary attractors. As this figure shows, the *a* to *b* transition has a lower MFPT compared to the other.

At each network state update *D*(*ŝ*) there is a probability that the state will change as a function of the Hamming distance (*h*) between the current state and the subsequent state *ŝ*_{t+1}← *D*(*η*(*ŝ*_{
t
}*, r*)). *MFPT* models uniform gene expression noise by considering probabilistic bit flips at every possible state of the network and deriving the distribution of passage times from analysis of the corresponding Markov process. Statistically, the probability distribution of bit flips can be seen as a binomial distribution, thus the probability of *r* bit flips, *η*(*ŝ*_{
a
}*, r*) is \left(\begin{array}{c}\hfill h\hfill \\ \hfill r\hfill \end{array}\right){p}^{r}{\left(1-p\right)}^{h-r}, where *p* is the probability of a single bit flip and *h* is the total number of bits.