Models of Cell Processes are Far from the Edge of Chaos

Complex living systems are thought to exist at the “edge of chaos” separating the ordered dynamics of robust function from the disordered dynamics of rapid environmental adaptation. Here, a deeper inspection of 72 experimentally supported discrete dynamical models of cell processes reveals previously unobserved order on long time scales, suggesting greater rigidity in these systems than was previously conjectured. We find that propagation of internal perturbations is transient in most cases, and that even when large perturbation cascades persist, their phenotypic effects are often minimal. Moreover, we find evidence that stochasticity and desynchronization can lead to increased recovery from regulatory perturbation cascades. Our analysis relies on new measures that quantify the tendency of perturbations to spread through a discrete dynamical system. Computing these measures was not feasible using current methodology; thus, we developed a multipurpose CUDA-based simulation tool, which we have made available as the open-source Python library cubewalkers. Based on novel measures and simulations, our results suggest that—contrary to current theory—cell processes are ordered and far from the edge of chaos.

Figure 1

Legend indicating model categories (marker shape) and specific highlighted models (marker color).

Figure 2

Comparison of four perturbation response measures (bold box borders) for a one-node oscillator. The unperturbed oscillator alternates between two states: its initial state $A$, which could be 0 or 1, and the opposite state, $\neg A$, which is 1 if the initial state is 0, and 0 if the initial state is 1. The perturbed trajectory begins with the oscillating node in the opposite state compared to the unperturbed trajectory, but otherwise its time evolution proceeds in the same fashion. At each time step $t$, the Hamming distance $h^t$ is computed. In the special case of $t=1$, $h^1$ is the Derrida coefficient $\delta$, which evaluates to 1 in this case. Indeed, $h^t=1$ for all $t$, so the asymptotic average of the Hamming distance, which we call the final Hamming distance (denoted $h^{\infty }$) evaluates to 1 as well. Alternatively, we can compute and compare the average behavior of the two trajectories. In both cases, the node is in the 0 state for half of the time steps, and in the 1 state in the other half. Thus, the average node value is 0.5 for both trajectories, and the fragility $\phi$, defined as the difference in these averages, is 0. Furthermore, we can consider a more coarse-grained averaging, where we compute the probability that a randomly perturbed node (in this example there is only one node to choose from) results in a different quasiattractor, i.e., a different pattern of fixed and oscillating nodes; the complement of this probability is a measure of robustness we call the quasicoherence. In this case, perturbing the initial state always results in the same quasiattractor (in which the sole node oscillates), so the quasicoherence is 1.

Figure 3

Distribution of update dependence in the Cell Collective. The root mean squared (RMS) difference between the node values when using synchronous or asynchronous update, as defined in Sec. 2a, is shown. The peak near zero indicates a high degree of timing robustness in the Cell Collective models. Representative models are indicated by symbols according to Fig. 1.

Figure 4

Short- and long-term perturbation responses in the Cell Collective measured in a phase-sensitive way. In the “Robust” regime (lower left quadrant) both short-term and long-term responses are below 1, which indicates perturbation extinction and is characteristic of ordered dynamics. In the “Sensitive” regime (upper right quadrant) both short-term and long-term responses are above 1. This indicates perturbation growth, which, in the extreme case, is characteristic of disordered or chaotic dynamics. The other two quadrants indicate cases of disagreement between the short-term and long-term responses. The short-term perturbation response $\delta$ has a slight correspondence with the long-term perturbation response under the specific setting when $h^{\infty }$ is monitored and synchronous update is used, in which the phase shifts are conserved. The relationship between short- and long-term responses is stronger when source nodes are fixed (right panel). The dashed line indicates the $y=x$ diagonal. The symbols indicate the model categories and highlighted models as defined in Fig. 1.

Figure 5

Scatterplot of the synchronous quasicoherences of the Cell Collective models when source nodes are ($x$ axis) or are not ($y$ axis) candidates for perturbation (the asynchronous distribution is available in Fig. 17 of Appendix pp6). When the values of source nodes are fixed, the quasicoherence values are tightly clustered around 1, indicating a high degree of phenotypic robustness. The symbols indicate the model categories and highlighted models as defined in Fig. 1.

Figure 6

Short- and long-term perturbation responses in the Cell Collective measured in a phase-insensitive way. In the “Robust” regime (lower left quadrant), both short-term and long-term responses are below 1, which indicates perturbation extinction and is characteristic of ordered dynamics. In the “Sensitive” regime (upper right quadrant), both short-term and long-term responses are above 1. This indicates perturbation growth, which, in the extreme case, is characteristic of disordered or chaotic dynamics. The other two quadrants indicate cases of disagreement between the short-term and long-term responses. In contrast with the traditional approach depicted in the left panel of Fig. 4, this figure illustrates perturbation response when source nodes and phase shifts are accounted for. Most models show a substantially more robust perturbation response when these factors are taken into consideration. The symbols indicate the model categories, and highlighted models as defined in Fig. 1.

Figure 7

Performance comparison of cubewalkers, cana, and booleannet on consumer hardware. 72 Cell Collective models were run using each tool using synchronous update. Timings were generated on a PC with an AMD Ryzen 5 3600X CPU at 3.8 GHz and a 2560 CUDA-core 1605 MHz NVIDIA 2070S GPU. Default methods were run without additional parallelization. For the cubewalkers tests, 2500 time steps and 2500 walkers (initial conditions) were used; for cana, 500 time steps and 500 walkers were used; and for booleannet, 100 time steps and 100 initial conditions were used. Thus, for each network, cana computed $5\times$ as many time steps for $5\times$ as many initial conditions as booleannet for an overall disadvantage of $25\times$. Similarly, cubewalkers computed $5\times$ as many time steps for $5\times$ as many initial conditions as cana, for a $25\times$ disadvantage relative to cana and a $625\times$ disadvantage relative to booleannet. The raw time to complete these tasks is plotted in the left panel, where we observe that cubewalkers consistently finishes its tasks an order of magnitude faster than the other methods, despite the fact that it has been given significantly more computational work. In the right panel, the average computation time per network node per time step per initial condition in these trials is plotted; this corresponds to the average (amortized) time to evaluate and apply an update function to a node. Here, we see that these amortized evaluations occur on the order of nanoseconds for cubewalkers, while they occur on the order of microseconds for cana and hundreds of microseconds for booleannet.

Figure 8

Performance comparison of cubewalkers, cana, and booleannet on a high-performance computer. Cell Collective models were run using each tool using synchronous update. Timings were generated using a workstation with two AMD EPYC 7542 CPUs (32 cores and 64 threads each) at 2.9 GHz and two 10 752 CUDA-core NVIDIA A6000 GPUs with 48GB of GDDR6 memory (only one GPU was used for the benchmarks). For the cubewalkers and cana tests, 2500 time steps and 2500 walkers (initial conditions) were used; for booleannet, 100 time steps and 100 initial conditions were used. For cana and booleannet, initial conditions were simulated in 128 parallel threads. On specialized hardware taking full advantage of parallelism, we see that the performance gap between cubewalkers and the other methods is narrowed compared to the performance gap on consumer hardware. Nevertheless, the gap remains considerable.

Figure 9

Standard deviation in the node average values as the number of walkers increases for the Cell Collective models. The three panels correspond to three different stages of the model simulation. The observed standard deviations agree well with the expectation based on Bernoulli random variables (continuous line). We chose the number of walkers such that the standard deviation is less than 0.01 (dashed vertical line).

Figure 10

Example networks and their state transition graphs (STGs) that illustrate update dependence. Panel A illustrates delay-dependency in the example of Eqs. (D1) and (D2), which are inspired by the Cell Cycle Transcription by Coupled CDK and Network Oscillators model [36]. Panel B demonstrates how much of the STG in the Metabolic Interactions in the Gut Microbiome model [40] is robust to changes in update scheme. Panels C–E illustrate how the asynchronous update can mix the synchronous attractor basins in a core regulatory circuit in the Colitis-associated colon cancer model [41] (panel C), the full synchronous STG of the Cortical Area Development model [42] (panel D), and a reduced version of the Apoptosis Network model [43] (panel E). In each interaction network, each node symbol contains the update function of the node. Blue edges ending in filled circles indicate positive regulation, and red edges ending in open circles denote negative regulation. In the STGs, attractor states are indicated by thick borders. The basin of attraction of each attractor is highlighted by the same color as the attractor. In asynchronous update, states can reach more than one attractor; such states are shaded using a gradient. In the lower part of panel A, states that differ only in the value of the delay node $Y$ are grouped together in shaded boxes.

Figure 11

Fragility of a four-node reduced version of the Human Gonadal Sex Determination model of [48]. Panel A depicts the interaction network. Each node symbol contains the update function of the node. Blue edges ending in filled circles indicate positive regulation, and red edges ending in open circles denote negative regulation. Panel B shows the state transition graph under the synchronous update. Attractor states are indicated by thick borders. The basin of attraction of each attractor is highlighted by the same color as the attractor. State transitions are shown with black arrows, and orange double-sided arrows indicate state pairs that are related by single-node perturbations. These are the transitions that can arise from single-node perturbations and that lead to different long-term behavior than is observed without perturbation. The thickness of each orange edge indicates the Hamming distance between the corresponding attractors. Panel C shows how to calculate the fragility of this reduced model exactly using the information in panel B.

Figure 12

The distribution of the Cell Collective models based on the number of source nodes (top) and the ratio of the number of source nodes to the total number of nodes (bottom).

Figure 13

Comparison of key measures for the 18 models in the Cell Collective that were altered to attain a better agreement with the originally published models.

Figure 14

Systematic evaluation of the dependence of the Derrida coefficient $\delta$ on the update scheme and on source node perturbations. The ensemble of Cell Collective models shows a general agreement between the Derrida coefficients obtained for synchronous and asynchronous update (top panels). When source nodes not candidates for perturbation, the Derrida coefficient dramatically decreases (bottom panels). For example, note that three cancer drug models (plus signs) lie far from the diagonal in the lower two panels, indicating that these models are highly affected by perturbations to source nodes. This is to be expected, as the source nodes in these models represent known cancer drugs that were selected because they have a tremendous impact on the behavior of cancer cells.

Figure 15

Relationships of the Derrida coefficient $\delta$ with the final Hamming distance $h^{\infty }$ and the fragility $\phi$.

Figure 16

Comparison of different ways to compute $h^{\infty }$. The ensemble of Cell Collective models shows an overall agreement between the final Hamming distances obtained for synchronous and asynchronous update (top panels). Exceptions include models that exhibit significant phase shifts under synchronous update. When source nodes are not candidates for perturbation, the final Hamming distance dramatically decreases (bottom panels).

Figure 17

Comparison of different ways to compute $q$. There is a general agreement between the quasicoherences obtained for synchronous and asynchronous update (top panels). Making the source nodes not candidates for perturbation dramatically decreases the fragility (bottom panels).

Figure 18

Comparison of different ways to compute $\phi$. There is a general agreement between the fragilities obtained for synchronous and asynchronous update (top panels). Making the source nodes not candidates for perturbation dramatically decreases the fragility (bottom panels).

Figure 19

Comparison of the final Hamming distance $h^{\infty }$ and the fragility $\phi$. The final Hamming distance is always larger than or equal to the fragility in both update schemes. The difference is much more prominent in the synchronous update in which any phase shift is permanent, compared to the asynchronous update in which the stochasticity can disperse it.