Approximating network dynamics: Some open problems

Discrete-time finite-state dynamical systems on networks are often conceived as tractable approximations to more detailed ODE-based models of natural systems. Here we review research on a class of such discrete models $N$ that approximate certain ODE models $M$ of mathematical neuroscience. In particular, we outline several open problems on the dynamics of the models $N$ themselves, as well as on structural features of ODE models $M$ that allow for the construction of discrete approximations $N$ whose predictions will be consistent with those of $M$.


1.
Introduction. Dynamics on networks can be complex. This is the main reason why there is so much interest in studying them. But it is also the main reason why so much of the literature on the subject relies on simulations. When it is possible to obtain analytic predictions, this is usually accomplished by approximating the dynamical system under study with a simpler, analytically tractable one, deriving predictions for the latter, and then showing that relevant features of the dynamics are not significantly altered by the approximation.
The purpose of this paper is to present some open problems that originate from prior results on such approximations by the author and his collaborators. In Section 2 we briefly review a result in mathematical neuroscience. This result inspired research on a certain class of finite-state dynamical systems on networks, and on more general questions of approximability of ODE systems by finite-state discretetime systems on networks. These topics will be discussed in Sections 3 and 4, respectively. The latter two sections are largely independent of each other. In both sections brief and far from exhaustive reviews of some existing results and open questions are given, with the goal of highlighting a few problems that in the author's opinion are likely to lead to some advances in the methodology for studying network dynamics.

2.
A motivating question: Dynamic clustering in neuronal networks. Empirical studies of firing patterns in some neuronal tissues have found that time seems to be partitioned into discrete episodes of roughly equal lengths with fairly sharp boundaries. For the duration of one such episode, a group of neurons fires together, while the other neurons are resting. Membership in the groups of firing or resting neurons may change from episode to episode, so that certain neurons sometimes fire together, and sometimes not. This puzzling phenomenon is known as dynamic clustering (see, for example, [23,26,30,31] and references therein).
Of course, discrete-time models of neuronal networks such as the ones described in Section 3 below would exhibit this kind of dynamics. But these models cannot provide any kind of explanation of the phenomenon of dynamic clustering, as the episodes are built right into the assumption of a discrete time line.
In [30] a class M of ODE models M for neuronal networks with a certain architecture was introduced. The architecture requires two layers of neurons: A layer E of excitatory neurons, and a layer I of inhibitory neurons. Each neuron is characterized by two variables, a voltage that evolves on a fast time scale, and a channel gating variable that evolves on a slow time scale. The ODEs for all neurons in E are identical, and so are the ODEs for all neurons in I. Each neuron in E takes synaptic input from a subset of I. There are no direct synaptic connections between any two neurons in E. On the other hand, each neuron in I takes synaptic input from a subset of E and all neurons in I.
In the models M ∈ M, time is continuous and not a priori partitioned into discrete episodes. No individual neuron acts as a designated pacemaker. Nevertheless, simulation studies reported in [30] show distinctive patterns of dynamic clustering for the neurons in E. Thus the architecture of these models, together with the intrinsic and synaptic properties of each neuron, may be considered a plausible explanation for the phenomenon of dynamic clustering.
One can now define a class N of discrete-time finite-state models N (see Definition 3.1 of Section 3), which could be considered approximations of the ODE models M ∈ M. The simulation results reported in [30] suggest that for each M ∈ M there exists some N ∈ N that will correctly predict, for all initial conditions in a large region U of the state space of M and all times T ≥ 0, which neurons in E will fire during which episodes. Clearly, in this case we would have a certain kind of consistency between the predictions of the models M and N . A more extensive discussion and formal definition of this notion will be given in Subsection 4.2.
The main theorem of [30] analytically confirms these simulation results. We give a simplified version of the theorem here; for details see Subsection 4.3 of [30].
Theorem 2.1. For each ODE model M ∈ M, if the intrinsic and synaptic properties of the cells are chosen appropriately, the dynamics of M will exhibit dynamic clustering. That is, time will be partitioned into episodes of approximately equal duration, with fairly sharp boundaries. During each episode, some group of neurons in E will fire, while the remaining neurons in this layer will be at rest.
Moreover, there exists a discrete-time finite-state model N ∈ N so that the models M and N are consistent.
As global dynamics of ODE models are notoriously difficult to study analytically, the simpler discrete-time models N ∈ N might be a welcome tool for deriving analytic predictions about the dynamics of the less tractable but more biologically realistic models M ∈ M. Provided, of course, that analytic predictions about the dynamics of N are feasible. As we will see in Section 3, the study of the models N already leads to a number of interesting mathematical questions.
A different set of questions naturally arises if we examine the nature of the consistency between the models M and N of Theorem 2.1. How should one define, in general, consistency between an ODE model M and a discrete-time approximation N ? And what structural features of ODE systems do, in general, favor the existence of discrete-time finite-state approximations as in Theorem 2.1? We will explore these questions in Section 4.

3.
A class of simple models for some neuronal networks.
3.1. Notation and definition. For the purpose of this paper, a network will be a directed graph, aka digraph, D with vertex set V D = [n] = {1, . . . , n} and arc set A D . A (discrete-time finite-state) dynamical system on the network is then defined by first specifying finite sets St i of states that node i can take at a given time t ∈ N = {0, 1, 2, . . . }. The state s of the system at time t is a vector s(t) = (s 1 (t), . . . , s n (t)), where s i (t) ∈ St i for all i ∈ [n]. Next we need to specify an updating rule for the change of s(t) over time. The updating rule must comply with the requirement that the state s i (t + 1) can depend only on the states s j (t) of those nodes j for which there is an arc j, i ∈ A D with source j and target i. In models with synchronous updating, the updating rule is applied to all variables simultaneously; in models with asynchronous updating, the variables are updated sequentially, in a fixed or random order.
Here we consider a class of models that mimic the firing patterns of certain networks of neurons that can induce firing in their postsynaptic neighbors, but need to go through a refractory period before firing again.
Definition 3.1. We let N be the class of discrete-time finite-state dynamical systems with synchronous updating that for every positive integer n, every loop-free digraph D = ([n], A D ) (that is, digraph with no arcs of the form i, i ), and every pair of vectors p = (p 1 , . . . , p n ), th = (th 1 , . . . , th n ) of positive integers contains exactly one neuronal network N on D with refractory periods p and firing thresholds th that is defined by St i = {0, 1, . . . , p i } and the following updating rules: (R) If s i (t) < p i , then s i (t + 1) = s i (t) + 1. (F) If s i (t) = p i and there exist at least th i distinct j ∈ [n] with s j (t) = 0 and j, i ∈ A D , then s i (t + 1) = 0. (E) If s i (t) = p i and there are fewer than th i nodes j ∈ [n] with s j (t) = 0 and j, i ∈ A D , then s i (t + 1) = p i .
In terms of neuroscience applications, s i (t) = 0 represents the event that neuron i fires at time t. Condition (R) formalizes a requirement that each neuron fires for exactly one time step and has to go through a refractory period of at least p i time steps before it can fire again. Once neuron i has reached the end of its refractory period, it will fire if it receives firing input from at least th i presynaptic neurons (condition (F)), and remain at the end of its refractory period otherwise (condition (E)).
The complete version of Theorem 2.1 above that appeared in [30] shows that for each N ∈ N there exists an ODE models M ∈ M such that M and N are consistent.
The dynamics of neuronal networks N ∈ N has been studied in a number of publications. An accessible introduction that reviews early results and contains projects suitable for REUs is given in [19]. But even for relatively simple fixed connectivities, such as a Hamiltonian cycle with one shortcut, characterization of all possible dynamics is far from easy [3,4]. As the complete connectivities of large actual neuronal networks are not currently known and may not be entirely hardwired by the genetic code, much attention has been given to studies where D is an Erdős-Rényi random digraph. In [18], two distinctive phase transitions in the dynamics of N ∈ N for relatively large mean indegrees were characterized. More precisely, when the mean indegree is above a threshold that scales like log(n) and the system starts in a typical initial state, all nodes will always fire as soon as they reach the end of their refractory periods. Right below this threshold, this will remain true for the vast majority of nodes, but no longer for all nodes. The latter will remain true as long as the mean indegree stays above a certain constant c > 1, which represents the threshold for the second phase transition. The exact value of c depends on the particular distribution of refractory periods and firing thresholds. The case of a third phase transition when the mean indegree is ≈ 1 was studied in [16,17]. The journal version [17] is restricted to the case when p = th = (1, 1, . . . , 1); its findings will be reviewed in Subsection 3.2. The preprint version [16] contains a number of generalizations of these results to arbitrary refractory periods p, and some results about arbitrary firing thresholds. Generally speaking, when th = (th, th, . . . , th), one would expect the third phase transition to occur when the mean indegree is ≈ th, but precise characterizations along the lines of the results reviewed in Subsection 3.2 remain open research problems at this time.
It is of interest to study the dynamics of N for other types of random connectivities, such as multitype Erdős-Rényi random digraphs, generic scale-free networks, and small-world networks. Research on these questions is currently ongoing [21]. It is also of interest to study the dynamics of models with updating rules that are similar to, but not identical with, the ones of Definition 3.1 and may be more biologically realistic for certain types of neuronal networks. Results for some modified updating rules have been reported in [5,12,29]. In the models of [5], both the excitatory neurons in E and the inhibitory neurons in I contain additional calcium ionic currents and there are more complex interactions among the cells. A reduction of the resulting ODE models M to discrete-time dynamical systems N is given in that paper. In N , a cell can fire for several consecutive episodes before other neurons take over. Simulations show dynamic clustering and suggest some form of consistency between M and N , but no rigorous result along the lines of Theorem 2.1 is currently known for these models. In [12] a more general class of networks than N of Definition 3.1 is studied. The models studied in [29] are very similar to the ones considered here, but nodes will fire for a fixed number b > 1 time steps before entering their refractory periods. [17]. The basic setup is as follows: We assume that the connectivity D of a neuronal network N ∈ N is randomly drawn from the distribution of Erdős-Rényi random digraphs with vertex set [n] so that each potential arc i, j is included in A D with a given probability π(n). The initial states s(0) are randomly drawn from the uniform distribution of all possible states in the resulting network N . For all i ∈ [n] we set p i = th i = 1.

Review of results from
Since the state space of N is finite and the updating rule is deterministic, each trajectory must eventually reach an attractor, that is, a set of states Att that it will visit infinitely often, in fixed order. Let α denote the length or size |Att| of this attractor, and let τ denote the length of the transient (part of the trajectory of s(0)), that is, the number of states that are visited only once. The focus of [17] is on deriving scaling laws for α and τ as n → ∞ for certain choices of π(n). These predictions will be presented below by expressions like: "τ = Θ(log n) a.a.s." Here a.a.s. abbreviates "asymptotically almost surely" and the statement means that there are constants C 2 > C 1 > 0 such that, with probability approaching 1 as n → ∞, the inequalities C 1 ln n < τ < C 2 ln n will hold. To assert only the first (second) of these inequalities, we write "τ = Ω(log n) a.a.s." ("τ = O(log n) a.a.s.").
Of particular interest is the question whether π(n) can be chosen in such a way that some form of complex dynamics becomes generic in the corresponding networks. The networks studied in [17] are examples of Boolean networks, and it makes sense to adopt the conceptualization of complex dynamics proposed by S. A. Kauffman and his followers [22] (see also [6,10] for more recent reviews), where they occur near the boundary between the ordered and the chaotic regimes. In the ordered regime, the lengths of attractors and transients tend to scale polynomially with network size, whereas in the chaotic regime they tend to scale exponentially. Thus scaling laws for α and τ that would simultaneously be superpolynomial and subexponential would provide strong evidence of the complex regime for the corresponding choices of the parameters π(n). Identifying regions of the parameter space where complex dynamics are predicted is relevant to neuroscience, as a number of results in the literature indicate that the potential of networks to perform sophisticated computations is highest if their dynamics falls into the complex regime (e.g., [14,28] and references therein). Another aspect of neuroscience applications is discussed in [4].
Simulations of the dynamics when π(n) = c n for some constant c suggest highly ordered dynamics when c < 1 or c > 1. They also suggest that the most complex dynamics occur when c ≈ 1. The main part of [17] backs up these empirical findings with rigorously derived scaling laws for α and τ that we will now briefly review. Some generalizations of these findings to neuronal networks with higher refractory periods and firing thresholds were reported in [16]. The next two results summarize Theorems 3.2 and 3.3 of [17], where Parts (c) and (d) of Theorem 3.3 are more precise versions of Theorem 3.4(a) of [17] that can be extracted from its proof.
Intuitively, one might expect a phase transition in the dynamics when π(n) ≈ 1 n , as this is when a giant strongly connected component that comprises a fixed fraction of all nodes appears a.a.s. in the underlying Erdős-Rényi digraph D. Theorem 3.2 shows that both above and below this phase transition the dynamics will fall into the ordered regime, with the lengths of attractors being bounded from above with high probability by a constant, and lengths of transients scaling a.a.s. at most polynomially. This confirms our empirical observations from simulations. Theorem 3.3 goes beyond these observations though in the sense that it clearly demonstrates that in a part of the lower end of the critical window where the phase transition in the network connectivity occurs, the lengths of the attractors will scale simultaneously superpolynomially and subexponentially. This can be taken as a clear indication of genuine complex dynamics in this region of the parameter space.
These results leave a number of open problems that are discussed in some detail in [17]. Here let us just focus on the two that in the author's opinion have the greatest potential for leading to methodological innovations. Note that Theorem 3.2 only tells us that 1 ≤ c crit ≤ 2, but simulations strongly suggest that c crit = 1. Similarly, Theorem 3.3 indicates that at the lower end of the critical window the lengths of transients will scale only sublinearly, which is not exactly what one would expect if the dynamics were truly complex.
for all k ≤ 2 , and the updating rule of the network dynamics implies that all nodes in C will be minimally cycling, that is, will fire at every other time step. It was then shown that, in a sense, these directed cycles C will entrain all other nodes in the giant strongly connected component of D so that they become minimally cycling from some time step on. Standard estimates on the diameter of the components of D then imply that this will a.a.s. happen relatively quickly, which gives the upper bound on the transients when c > 2.
Notice that this argument relies on the presence of certain small-scale structures in a randomly chosen initial state. In order to make progress on Problems 1 and 2 though, it seems necessary to explicitly consider distributions of states at some later time. In particular, for a positive answer to Problem 1 one would want to show that a.a.s. there will exist a directed cycle C as in the previous paragraph such that at some time t > 0 we have s i k+1 (t) = 1 − s i k (t) for all k ≤ 2 . This would require some knowledge of the distributions of the structures (A D , s(t)) for time t > 0.
Let us look at this approach from an abstract perspective that applies to all discrete-time dynamical systems on networks. The distributions of initial states s(0) and network connectivities A D uniquely determines, for every t > 0, a distribution of s(t), which can be treated as a random coloring of the nodes of a given random digraph (V D , A D ). These distributions of the random colorings s(t) will change over time, in a way that is determined by the connectivity A D and the updating rule. It seems intuitively clear that knowledge about the distributions of s(t) (or of the pairs (A D , s(t)) if the connectivity is also randomly drawn) might allow us to derive conclusions about important features of the network dynamics. Therefore the author of this note believes that methods for studying the distributions of (A D , s(t)) might be very useful for studying a wide variety of networks dynamics. Alas, these distributions seem to be as difficult to study as the network dynamics itself.
However, the situation is not entirely hopeless and at least some meaningful inferences about these distributions are feasible. For a node i of the network with connectivity D, define its (backward) neighborhood of depth as the set N(i, ) of all nodes j such that there exists a directed path from j to i of length ≤ , and let D(i, ) be the subdigraph of D that is induced by N(i, ). Now note that for all dynamical systems as described at the beginning of this section, the distribution of the sequencess i ( ) = (s i (0), s i (1), . . . , s i ( )) is determined entirely by the distributions of the subdigraphs D(i, ) and the restrictions of initial states to N(i, ). In many cases the latter distributions may have relatively simple approximations; we will give some examples shortly. Moreover, if the distance between nodes i and j is larger than , then for anyā,b the events thats i ( ) =ā ands j ( ) =b are independent. Thus the overall distribution of the sequencess i ( ) can be approximated using the Chen-Stein method (see [7] for a book-length treatment). This method has already been used in the literature on network dynamics [13], but only sporadically.
But back to our neuronal networks. For the proof of Theorem 3.3(c)(d), we analyzed the structure of Erdős-Rényi random digraphs in the lower part of the critical window that is covered by the assumptions. For these connection probabilities the network will have a.a.s. many small strongly connected components that contain one directed cycle each, of lengths q r r , where r are distinct prime numbers. In a randomly chosen initial condition, a.a.s. the vast majority of these directed cycles will contain at least one node in the resting state 1 and at least one node in state 0, which is interpreted as firing. We then showed, essentially by considering the dynamics on D(i, ) for i in such a cycle and sufficiently large, that a.a.s. the firings in most of these directed cycles will perpetuate indefinitely. This observation was then used to derive a lower bound of r r on the length of the attractor. This argument was made possible by the structure of the connected components at the lower end of the critical window, which are known to a.a.s. contain at most one directed cycle. For progress on Problem 2, one would like to extend this approach to models where, for example, π(n) = 1 n and try to derive lower bounds on the lengths of the transients. However, for this one would need to consider induced subgraphs D(i, ) that contain multiple directed cycles.
A similar approach may lead to progress on Problem 1. For a given π(n) ≤ c n and sufficiently large n, a.a.s. for the vast majority of nodes i the subdigraph D(i, ) will be a directed tree with root i. Moreover, for 0 < t ≤ , we will have s i (t) = 0 if, and only if, there is a node j ∈ N(i, ) such that the unique directed path (j = j 0 , j 1 , . . . , j t = i) from j to i in D(i, ) has length t and the property that s jtt (tt) = 0 for all 0 ≤ tt ≤ t. If this is the case, then we say that the firing of node i at time t is induced by node j. This observation shows that for fixed t, large n and the vast majority of nodes, the distributions of firing patterns (s i (0), s i (1), . . . , s i (t)) will be given by the analogous distributions for the root of a random tree that can be generated by a Galton-Watson birth process with mean number of offspring π(n)(n− 1). Working out the latter distributions should be feasible. Feasible, most likely yes, but not entirely easy. At the time of this writing, the author has only obtained some partial results on these distributions. The interesting situation occurs of course if the Galton-Watson process that generates D(i, ) does not die out prematurely and generates a tree of height + 1. When π(n) ≥ 1 n , this will be the case for the nodes i in the largest strongly connected component of D for which D(i, ) is indeed a tree. It is essentially already shown in one argument of [17] that then for all t − > 0 and all sufficiently large t and n, a.a.s. for the vast majority of nodes i in this largest component there will be some stretch of the form 1, 0, 1, . . . , 0, 1).
The time t = t(t − ) required for appearance of such a stretch of uninterrupted firings of node i to appear will depend on π(n). Now we might consider a related digraph D − = D − (t 0 , t − ) whose nodes are restricted to those for which we observe a stretch of uninterrupted firings from time t 0 until time t = t 0 + t − , and whose arcs are of the form j, i , where the firing of node i at time t is induced by the firing of j at time t 0 . Thus at time t 0 any directed cycle in this digraph would show the same firing pattern as the directed cycle C in our argument for c > 2 that was sketched above. Thus if we could prove that D − (t 0 , t − ) contains a directed cycle (necessarily of even length), then we would obtain a positive answer to Problem 1.
In general, D − (t 0 , t − ) will not be an Erdős-Rényi digraph, as the events j, i ∈ A D − and j , i ∈ A D − may be dependent. However, if i , j / ∈ N(i, t) ∪ N(j, t) or i, j / ∈ N(i , t) ∪ N(j , t), then these events will be independent. Thus it may be possible to derive global properties of D − (t 0 , t − ) by analyzing the subdigraphs D(i, 2t) or by using the Chen-Stein approximation. At the time of this writing the author does not know what fraction of nodes in the largest strongly connected component would be expected to be nodes of D − (t 0 , t − ), but conjectures that when π(n) = c n for some c > 1 and t is large enough, this fraction will be close to 1.

4.1.
More motivating examples: Models of gene regulation. The set of biochemical reactions that will occur in a cell at a given time depends on the enzymes that are available for catalyzing these reactions. These enzymes in turn are proteins coded by certain genes. Thus the biochemical state of a cell is determined by which genes are being expressed at the time. The possible expression patterns form the states of gene regulatory networks.
The dynamics of gene regulatory networks can be conceptualized in a variety of mathematical frameworks (see, for example [8,27]), including ODE models. Since it is difficult to measure concentrations of the relevant biomolecules with reasonable accuracy, much work has been done on modeling gene regulatory networks with Boolean networks, where the concentration is discretized into either Boolean value 0 (low) or Boolean value 1 (high). While ODE models appear to be biologically more realistic, in many cases it has been found that Boolean models make fairly accurate predictions and can elucidate relevant mechanisms of gene regulation. A classical example is the Boolean model of the yeast cell cycle network of [24].
There may be important biological reasons why Boolean models work as well in this context as they apparently do. These models are based on the idea of binary switches, and cells may need to implement reliable switches for proper functioning. For example, the cell cycle involves several distinctive stages, with characteristic gene expression patterns during each stage, and fairly rapid switches between these patterns when the cell moves to the next stage. For proper functioning in the face of the stochasticity inherent in biochemical reactions, the underlying continuous-time dynamics must be sufficiently robust to exhibit these switches for initial conditions in large regions of the state space. This robustness may make it feasible to study gene regulatory networks in the Boolean formalism.
These observations suggest that reasonably realistic ODE models M of gene regulation might be consistent with Boolean approximations N . Indeed, a number of results along these lines have been reported in the literature. For example, Leon Glass and his followers have produced a large body of work on consistency (with asynchronous, but not synchronous updating) between piecewise linear ODE models of gene regulatory networks and their Boolean approximations; see [11] for a review. Gehrmann and Drossel [15] proved consistency (with synchronous updating) between an ODE model M and a Boolean model N for a certain small gene regulatory network. Simulation studies reported in [1] indicate that these results can be extended to larger gene regulatory networks. However, inconsistencies between ODE models of gene regulatory networks and their discrete-time approximations have also been reported; see, for example [25] and references in [27].
Interestingly enough, in terms of the biochemistry of the cell, what appears to really matter is reliable switching between concentration levels of gene products (such as enzymes or transcription factors whose binding to certain regions of the DNA can initiate or terminate expression of a gene). However, gene expression is a two step process of first transcribing the code for a gene into messenger RNA (mRNA) and then translating the mRNA into the gene products. Thus in gene regulatory networks, gene products play a role somewhat similar to the neurons in layer E of the neuronal networks of Theorem 2.1, while the mRNA mediate their interactions and play a role somewhat similar to the neurons in layer I of these neuronal networks. These observations suggest that variables that play the role of intermediaries may be beneficial or even necessary for consistency between certain ODE models and their Boolean approximations. Moreover, the model of [24] and similar ones suggest that robustness of switching behavior is achieved by multiple feedback loops in the network. The effect of such feedback loops has been studied from a mathematical perspective in a number of publications; see, for example, [2,9,27] and references therein. However, more realistic ODE models than the piecewise linear ones of Glass and his collaborators are difficult to analyze, and general conditions for consistency between these models and their Boolean approximations appear difficult to derive.

4.2.
Consistency between flows and Boolean systems. Consider a system of autonomous ODE's with variables x = (x 1 , . . . , x m ), so that for each i ∈ [m]: Time can take any real values here and is denoted by T to distinguish it from time t in discrete models where it takes values in N. Under suitable assumptions about the functions F i , (1) will generate a differentiable flow M on a compact manifold K, with trajectories defined at all times T ∈ R. For convenience, we will refer to both the system (1) and the corresponding flow as an ODE model M .
We can now define a digraph D on E by including j, i in the arc set A D if, and only if, there exist Boolean states s, s that differ only in variable j (so that s j = s j and s k = s k for all k = j) with f i ( s) = f i ( s ).
We fix a discretization S that maps real vectors x to Boolean vectors s = S( x) ∈ {0, 1} E . Typically, this function will be defined in terms of thresholds Θ i for i ∈ E so that S( x) i = 1 if, and only if, x i ≥ Θ i . Any such discretization assigns to each forward trajectory x(·) : [0, ∞) → K a function S( x(·)) : [0, ∞) → {0, 1} E , defined by S( x(·))(T ) = S( x(T )). We will call it a real-time Boolean trajectory. (b1) T 0 = 0 and T t < T t+1 for all t ∈ N.
(c) If (b5) can be strengthened, for all x(0) ∈ U , to (b5) If s(t + 1) = s(t), then s(t + 1) = f ( s(t)), then we say that the models M and N are strongly consistent.
Note that consistency in the sense of Definition 4.1 corresponds to an intuitive notion of "consistency with asynchronous updating," while strong consistency corresponds to "consistency with synchronous updating." Theorem 2.1 asserts strong consistency when the times T t are taken as the midpoints of the "episodes" that were mentioned at the beginning of Section 2. The theorem works both for the Boolean case that is explicitly covered by Definition 4.1 and its straightforward generalization to models N whose variables can take more than two states.
There are obvious similarities between the notions of consistency that we are interested in here and digital computations. Every implementation of a numerical ODE solver implicitly approximates real-time real-valued trajectories of an ODE model with initial segments of discrete-time trajectories in a finite-state model. Much research in numerical analysis focuses in essence on questions of how reliable these approximations are over a long but finite time horizon for given ODEs that may be somewhat ill-behaved. In contrast, (strong) consistency requires that certain flows M are sufficiently well-behaved so that the discrete-time system N makes correct predictions about the trajectories of the flow M at all positive times. The author of this paper believes that it is of interest to study, at a general level, which structural features of ODE models M favor or imply such consistency or strong consistency.

4.3.
A class of toy models. In order to gain insight into general mechanisms that favor consistency or strong consistency with Boolean models, two classes D 1 and D 2 of toy ODE models were introduced and studied in [20]. The dynamics of the flows M for these systems are relatively easy to understand, and each class is rich enough so that every Boolean network N with n variables becomes a natural candidate for a Boolean approximation to one of these flows.
More precisely, let N be a Boolean network with updating functions f i : {0, 1} [n] → {0, 1}. We then convert N into an ODE model M according to a fixed conversion scheme. These constructions are partially based on the ones proposed in [32]. In the interest of greater transparency, here we only outline the basic mechanism; for the technical details we refer the reader to [20].
When N is converted into M ∈ D 1 , the variables x i in M are enumerated by the same set [n], so that each of them has a Boolean counterpart. A Boolean discretization S of the state space of M is chosen as described in Subsection 4.2, with threshold Θ i = 0 for all i. The important feature of our models M is that if we fix the values of all parameters x j for j = i, then we can relatively easily study the intrinsic dynamics M i of variable i. This corresponds to reducing the ODE system of M to the equation for dxi dT that is actually shown in (1) and treating all other values x j as fixed parameters. The combined effect of these parameters is determined by a real value Q f i ( x) that depends only on those x j for which j, i is an arc in the connectivity of N .
, then there exists a globally stable equilibrium x * * i > 0 in M i ; and if Q f i ( x) is moderate, then in M i there exist two locally stable equilibria x * i < 0 < x * * i together with an unstable equilibrium. Thus Q f i ( x(t)) acts as a bifurcation parameter that changes over time, with two bifurcation values q * < q * * . The function Q f i embodies the above mentioned conversion scheme of f i . Essentially, in a state where all variables x j of M have sufficiently large absolute values, Q f i ( x) < q * when f i (S( x)) = 0 and Q f i ( x) > q * * when f i (S( x)) = 0.
The construction of models M in D 2 is similar, except that each model M ∈ D 2 has m = 2n variables. We call the variables x i with i ∈ [n] signature variables; they have Boolean counterparts in N . The remaining variables in M are called signaling variables. The dynamics of x i+n in M is governed by the same ODE as for the variable x i in a corresponding model in D 2 . In particular, it is directly influenced only by the signature variables. In contrast, the ODE for each signature variable i depends only on x i+n . Thus interactions between the signature variables are mediated by the signaling variables. The latter play a role somewhat similar to mRNA concentrations in gene regulatory networks and evolve at slower time scales than the signature variables.
In order to investigate consistency between M and N for these constructions, we considered the set of switching times for x(0), defined as follows: Here SW ( x(0)) is always a closed subset of [0, ∞). We can then define a sequence (T t ) t∈N = (T t ( x(0))) t∈N by taking T j as the midpoints of the intervals between two successive switching times if there is no last such time, or in terms of fixed increments after such a last switching time. In [20], M and N were considered (strongly) consistent if the resulting sequence satisfied the requirements of Definition 4.1.
This conceptualization of the sequences (T t ) t∈N implies that for generic trajectories the Boolean value S( x i (t)) of at most one signature variable can change in a given interval (T t , T t+1 ). Thus this framework permits strong consistency only if the updating function f = (f 1 , . . . , f n ) of N is what we called one-stepping, which means that for every s the Boolean vectors s and f ( s) differ in at most one coordinate. Hence the following theorem of [20] is in a sense the strongest possible one that can be obtained about strong consistency: Theorem 4.2. Let N be a Boolean system with a one-stepping updating function f and let M ∈ D 2 be a conversion of N . Then as long as the separation of time scales for the evolution of the signature variables and signaling variables in M is sufficiently large, the models M and N will be strongly consistent.
However, more general results can be obtained for consistency instead of strong consistency. For example, we also have a more technical notion of monotonestepping Boolean functions. All one-stepping Boolean functions are monotonestepping, but not vice versa. For precise details, see Definition 4 of [20]. Theorem 4.3. Let N be a Boolean system with a monotone-stepping updating function f and let M ∈ D 2 be a conversion of N . Then as long as the separation of time scales for the evolution of the signature variables and signaling variables in M is sufficiently large, the models M and N will be consistent.
These results conform to expectations that one might form based on simulation studies for gene regulatory networks reported in [1]. However, the assumptions in these theorems appear to be overly restrictive. Indeed, we know that some additional assumption on f is needed in the Theorem 4.3, but being monotonestepping is not a necessary one. With regard to Theorem 4.2 we know that it is the strongest possible one for the definition of (T t ( x(0))) t∈N based on SW ( x(0)), but there might be a better way to conceptualize strong consistency.
Problem 4. Does Theorem 4.2 generalize to a larger class of Boolean networks N if the sequences (T t ( x(0))) t∈N are chosen by a different method than in [20]?
While the proofs of Theorems 4.2 and 4.3 given in [20] rely on the particular details of the intrinsic dynamics of the signature variables, the author of this note believes that they should remain valid for broader classes of ODE models. Recall the above partial description of the classes D 1 and D 2 in terms of bifurcations of the internal dynamics M i for the signature variables and their interactions in the global dynamics. The deliberate vagueness of this description is meant to suggest that all that may really matter in Theorems 4.2 and 4.3 is this general structure of bifurcations, where locally stable attractors that are contained within regions that get discretized to a fixed Boolean value may disappear, appear, or become globally stable in response to parameter changes that are mediated by the signaling variables that change on a slower time scale. This makes it highly plausible that some general theorems along the lines of Problem 5 can be proved. However, we want to caution against overinterpreting Problem 5 as a quest for finding the universal mechanism for (strong) consistency. In particular, the models in the class M of Theorem 2.1 are very different from the ODE models in D 1 and D 2 .
Interestingly enough, the analogues of Theorems 4.2 and 4.3 do not hold for the class D 1 ; see [20] for some examples. Thus signaling variables of sorts seem to be needed to ensure robustness (conceptualized here as consistency) of switching behavior, even when the Boolean updating function is one-stepping.
Could we get strong consistency in the sense of Definition 4.1 for arbitrary Boolean systems? Not with ODE systems M ∈ D 2 for the standard choice of signature variables. But an unpublished result from the work on [20] shows that every Boolean system N − can be embedded into a one-stepping Boolean system N with additional variables. By Theorem 4.2, there is an ODE model M that is consistent with N . If we treat only those variables x i of M that have a Boolean counterpart in N − as signature variables (the set E of Definition 4.1) and the remaining variables as signaling variables, then we do get strong consistency between M and N − .
Admittedly, this construction of M from N − is rather contrived and would require a very large number of signaling variables. Biological networks appear to achieve robustness of switching behavior of the relevant "signature variables," conceptualized here as (strong) consistency, by incorporating a variety of feedback loops. The papers [2,9,24,27] that were cited in Subsection 4.1 constitute only a small sample of the mathematical studies of this kind of network structure. Evolution, a stochastic process, needs to somehow discover architectures that will exhibit the required robustness, perhaps by randomly making and breaking connections in the network connectivity. Therefore one might expect that some version of consistency becomes a generic property when these feedback loops of the signaling variables are drawn from a suitable distribution. This brings us round circle to dynamics on networks with random connectivities and to our last problem: Problem 6. Explore notions of consistency for versions of the class D 2 that incorporates random feedback loops.