General method to find the attractors of discrete dynamic models of biological systems

Analyzing the long-term behaviors (attractors) of dynamic models of biological networks can provide valuable insight. We propose a general method that can find the attractors of multilevel discrete dynamical systems by extending a method that finds the attractors of a Boolean network model. The previous method is based on finding stable motifs, subgraphs whose nodes’ states can stabilize on their own. We extend the framework from binary states to any finite discrete levels by creating a virtual node for each level of a multilevel node, and describing each virtual node with a quasi-Boolean function. We then create an expanded representation of the multilevel network, find multilevel stable motifs and oscillating motifs, and identify attractors by successive network reduction. In this way, we find both fixed point attractors and complex attractors. We implemented an algorithm, which we test and validate on representative synthetic networks and on published multilevel models of biological networks. Despite its primary motivation to analyze biological networks, our motif-based method is general and can be applied to any finite discrete dynamical system.

Figure 1

Demonstration of the construction of a quasi-Boolean regulatory function. A 3-level node $A$ has regulatory function $f_A=B+C$, where $B$ and $C$ both have 2 levels. From the truth table, one can identify the regulatory function for each virtual node of $A$, by connecting all conjunctive clauses that yield the same state of $A$ with the Boolean “or” operator. In this way, each virtual node's regulatory function will have a Boolean disjunctive normal form.

Figure 2

Example of the multilevel Quine-McCluskey algorithm. A Boolean node $D$ is regulated by a Boolean node $A$ and two 3-state nodes $B$ and $C$. The original function of $D$ is shown in a truth table on the top left, in a form summarizing all input combinations that yield $f_D^{\left(1\right)}=1$. The top right table shows the minterms sorted according to the number of zeros in them. From this table, one can merge the terms between layers that are different by 1 digit, if all states of the difference node are present within the two layers. The result of the merging is shown below. Merged terms are represented by an $X$. There are 5 leftover terms after first-order merging, and there is 1 leftover term after zeroth-order merging. The sum of all six terms is the final expression.

Figure 3

Construction of an expanded network from a regulatory function. Virtual node $A0$ has function $f_A^{\left(0\right)}=B0$ or ($C1$ and $B1$), so in the expanded network, $B0$ is connected directly to $A0$; $C1$ and $B1$ are connected indirectly to $A0$ via composite node “$C1$ and $B1$”. $A1$ has function $f_A^{\left(1\right)}=C0$ and $B1$, so $C0$ and $B1$ are connected indirectly to $A0$ via composite node “$C0$ and $B1$”.

Figure 4

Illustration of stable motif identification in a three-node network. (a) The original network and the regulatory functions of each node. (b) The expanded network is constructed according to the steps in Sec. 1 E, and then the stable motifs are found by their definition in Sec. 1 F. (c) Stable motifs found in this example. The first stable motif, $A0$, $B0$, corresponds to a fixed point attractor of the system $A=0, B=0, C=0$. The state $C=0$ is found by plugging $A=B=0$ into the regulatory function of $C$. The second stable motif corresponds to another fixed point attractor $A=2, B=2, C=1$.

Figure 5

An example of an oscillating motif in a multilevel network. Panel (a) shows the network and regulatory functions; panel (b) indicates the expanded network and motifs. $A0$ and $B0$ form a stable motif, indicating a fixed point $A=0, B=0$, while $A1$, $A2$, $B1$, and $B2$ form an oscillating motif, indicating a possible complex attractor involving states $A=1, A=2, B=1$, and $B=2$. Panel (c) indicates the state transition graph of the system when using general asynchronous update. The stable motif and oscillating motif identified in (b) correspond to a fixed point and a complex attractor, respectively.

Figure 6

An example of an oscillating motif that contains a stabilized node. (a) The network and regulatory functions. (b) The expanded network and motifs. The oscillating motif contains only one virtual node of $B$, meaning that $B$ will stabilize at 1 in the complex attractor. (c) The state transition graph using general asynchronous update. There are two attractors: a fixed point attractor, and a complex attractor.

Figure 7

Attractor identification for a four-node network by a motif succession diagram. (a) The network and the regulatory function of each node. (b) Motif succession diagram. Three motifs are found from the original network, including 2 stable motifs ($A0$, $B0$), ($C1$, $D1$), and one oscillating motif ($A1$, $A2$, $B1$, $B2$). For each motif, the values of the nodes in the motif are plugged into the regulatory functions, reducing the network. Then new motifs are identified from the reduced networks. The sequences corresponding to the three motifs are labeled (1), (2), and (3).

Figure 8

Histogram of QM transformation runtime on 100 randomly generated heterogeneous networks with 50 nodes. The result shows that the complexity of QM transformation is much less than identifying motifs.

Figure 9

An example of a timing-dependent complex attractor. (a) The network and regulatory functions. (b) The state transition graph under synchronous update. Each node of the state transition graph is a state, given in the order $A, B$, and each edge is a state transition allowed by synchronous update. The system has two fixed points, (0, 0) and (1, 1). It also has a complex attractor formed by the states (0, 1) and (1, 0). (c) The state transition graph under general asynchronous update (i.e., when one node is updated at a time). Only the two fixed point attractors exist. The synchronous complex attractor is timing-dependent and does not exist in this update scheme.

Figure 10

An example of an oscillating motif without a complex attractor. (a) The network and regulatory functions. (b) The expanded network and motifs. There is a stable motif formed by $A0$ and $B0$, and an oscillating motif made up by $A1$, $A2$, $B1$. (c) The state transition graph using general asynchronous update. There is only one attractor, which is a fixed point. The transient oscillation between states (2, 1) and (1, 1) will eventually converge into the fixed point.

Figure 11

Examples of stabilized nodes downstream of oscillating node(s). (a) A Boolean example where $A$ and $B$ oscillate but their downstream $C$ is stable under that oscillation. (b) The general asynchronous state transition graph of nodes $A$ and $B$. The state ($A=1,B=1$) is not visited in the long term, leading to the stabilization of $C=0$. (c) A multilevel example where $A$ is oscillating between 1 and 2, leading to $B$ stabilizing at 1. This example arises because of asymmetry in the nodes’ number of states: $A$ has three states but $B$ only has two states.

Figure 12

An example of an unstable oscillation. The system has a fixed point and a complex attractor. (a) The network and regulatory functions. (b) The expanded network and motifs. The entire expanded network forms an oscillating motif, containing the stable motif by two nodes $A1$, $B1$, and one composite node. (c) The state transition graph using general synchronous update. There is a fixed point attractor $A=1, B=1$, and a complex attractor. Note that in the complex attractor, although both $A$ and $B$ are allowed to enter state 1, they cannot be in state 1 simultaneously.

Figure 13

An example of an oscillating motif containing two complex attractors. (a) The network and regulatory functions. (b) The expanded network and motifs. The entire expanded network forms an oscillating motif. (c) The state transition graph. For simplicity self-loops representing self-transitions are not shown in the graph. There are two complex attractors; the first attractor is $B=0, A=0$ or 1, and the second attractor is $B=1, A=2$ or 3.