Multi-parameter exploration of dynamics of 
regulatory networks
Author: Tomas Gedeon
© This manuscript version is made available under the CC-BY-NC-ND 4.0 license https://
creativecommons.org/licenses/by-nc-nd/4.0/
Gedeon, Tomáš. “Multi-Parameter Exploration of Dynamics of Regulatory Networks.” 
Biosystems 190 (April 2020): 104113. doi:10.1016/j.biosystems.2020.104113. 
Made available through Montana State University’s ScholarWorks 
scholarworks.montana.edu 
Multi-parameter exploration of dynamics of regulatory networks⋆
Tomáš Gedeona
aDepartment of Mathematical Sciences
Montana State University, Bozeman, MT 59715
ART ICLE INFO ABSTRACT
Keywords: Over the last twenty years advances in systems biology have changed our views on microbial
Boolean models communities and promise to revolutionize treatment of human diseases. In almost all scien-
network dynamics tific breakthroughs since time of Newton, mathematical modeling has played a prominent role.
parameters Regulatory networks emerged as preferred descriptors of how abundances of molecular species
depend on each other. However, the central question on how cellular phenotypes emerge from
dynamics of these network remains elusive. The principal reason is that differential equation
models in the field of biology (while so successful in areas of physics and physical chemistry),
do not arise from first principles, and these models suffer from lack of proper parameterization.
In response to these challenges, discrete time models based on Boolean networks have been
developed.
In this review, we discuss an emerging modeling paradigm that combines ideas from differ-
ential equations and Boolean models, and has been developed independently within dynamical
systems and computer science communities. The result is an approach that can associate a range
of potential dynamical behaviors to a network, arrange the descriptors of the dynamics in a
searchable database, and allows for multi-parameter exploration of the dynamics akin to bifur-
cation theory. Since this approach is computationally accessible for moderately sized networks,
it allows, perhaps for the first time, to rationally compare different network topologies based on
their dynamics.
1. Introduction
"Nowadays, many computational biologists avoid modeling-as-data-fitting, opting instead to create models in which networks are specified in terms of elements
and interactions (the network “topology”), but the numerical values that quantify those interactions (the parameters) are deliberately varied over wide ranges. As a
result, the study of such networks focuses not on the exact values of outputs, but rather on qualitative behavior, e.g., whether the network acts as a “switch,” “filter,”
“oscillator,” “dynamic range adjuster,” “producer of stripes,” etc. By investigating how such behaviors change for different parameter sets— an exercise referred to as
“exploring the parameter space”—one starts to assemble a comprehensive picture of all the kinds of behaviors a network can produce." (Lander 2004 [47])
Differential equations are inseparably linked to the birth of modern science. Among Newton’s fundamental con-
tributions to science was the introduction of a phase space where only those descriptors of a physical system that
have effect on its dynamics are represented. Second, the development of a concept of a differential equation that al-
gebraically links infinitesimal changes and other quantities allowed a quantitative prediction of a future state in the
phase space. Most important was his realization that, in the abstract world of the phase space, one can postulate the
existence of non-measurable quantities (i.e gravitational forces), solve a differential equation for measurable quantities
(i.e distances) and based on their comparison with experimental observations one is able to reject a particular form
of non-measurable quantity. This laid the foundation to the fruitful cycle of theory and experiments that has driven
science forward over the last 300 years.
Differential equations underwent a paradigm shift startingwith thework of Poincare at the start of 20th century. The
emphasis shifted from solving differential equations to qualitative theory, from a finite time horizon to the long term
behavior. The resulting tools from the qualitative theory of differential equations were instrumental to understanding
the wide range of phenomena from the singularities in n-body dynamics [74, 58] to the description of neural spike
propagation modeled by Hodgkin-Huxley equations [41].
⋆This research was partially supported by NSF TRIPODS+X grant DMS-1839299, DARPA FA8750-17-C-0054 and NIH 5R01GM126555-01.
In addition, part of this work has been done while the author was a visiting scientist at INRIA in Grenoble, France, whose support and hospitality
is graciously acknowledged. .
∗Corresponding author
gedeon@math.montana.edu (T. Gedeon)
ORCID(s): 0000-0001-5555-6741 (T. Gedeon)
T Gedeon: Preprint submitted to Elsevier Page 1 of 16
Dynamics of networks
Naturally, when the world of molecular and cellular biology opened up to experimental measurements, the toolbox
of differential equations had been readily available to try to model the dynamics of the intricate systems that underlie
cellular life. The problems that were encountered were both severe and fundamental; there seemed to be no easy fix
that could simultaneously address all of them.
In this paper we describe these difficulties and how a combination of Boolean networks and continuous time dy-
namics begin to converge to a modeling framework that is well suited for modeling regulatory networks that arise in
molecular and cell biology.
There are at least three fundamental challenges when using differential equations to model the dynamics of regu-
latory networks: nonlinearity, parameterization and continuity of phase space, parameter space and time.
1.1. Nonlinearity
Because of the multi-scale character of the biological networks, there is a lack of "first principles" description of
the nonlinear right hand side of the differential equations models. In order to analyze a differential equation, either
numerically or analytically, an explicit nonlinearity is necessary. Nonlinearities arise in physics from fundamental
laws like Newton’s law. There are no such well established laws in molecular biology that would govern the interaction
between genes, or between a protein and its receptor. What is being used instead is a parameterized nonlinearity in a
form of a Hill function, or nonlinearities that arise from mass action chemical kinetics.
1.2. Parameterization
Parameters pose another big problem. In Newtonian physics, the parameters are masses of objects that interact;
these are directly measurable. In contrast, in molecular biology the constants such as the "half-saturation rate", or a
“reaction rate", can only be measured by fitting a solution of a differential equation with an assumed nonlinearity to
the experimental data. If we assume a different nonlinearity, the value of the fitted parameter will be different. Such a
parameter can only be interpreted within the framework of its model.
This makes it difficult to permanently associate a parameter value to a given system. Therefore it is challenging
to build on the work of others since their parameters may have been used in a context of a slightly different model.
At best, one has to re-fit the data in the context of the new model; but more often than not, one must measure the
parameters again. Since measuring parameters in biology is expensive, parameter estimates are often poor. This leads
to parameterization of models by fitting the models solutions to the data. This results in identifiability problems, where
there are many different parameters and models that can fit the same data.
However, the exploration of the parameter space is often how the model is expected to inform biology: the interest
is to know what perturbation can bring the system from one state to another, what medical interventions, modeled as
a change of the parameters, have medically desirable effects, and how does underlying genetic diversity determine the
medical outcome. These are all questions about how the dynamics of the system changes with parameters.
1.3. Continuous time, phase space and parameter space
Continuous time, phase space and parameter space are arguably more realistic in modeling the natural world than
discrete time models. However, a continuous description is not directly representable in computers which always
brings up the issue of a proper approximation of a continuous system.
On the other hand, discrete models’ future behavior depends on the present state via a map, phase space is often
discrete (i.e. when we represent a state of a gene by being On or Off, as in the Boolean network models), or there
is only a discrete set of parameters. In a discrete description enumeration of all possibilities is, at least in principle,
feasible, and thus no approximation is necessary.
2. Models of gene regulatory networks
Even though without nonlinearities and parameters, the differential equation models cannot be expected to be very
useful, the demand for informative models in systems biology is growing. The potential medical applications demand
an analysis of large systems with dozens of interacting genes and proteins. These systems are often characterized in
terms of a network of interactions, which represent a direct up-or down-regulation of one gene, or protein, by another.
As a result, there is a great interest in understanding the dynamics of networks. However, the term itself provides
an insight to the conflict between expectations placed on mathematical models and its potential to deliver on them.
Biologists would like to understand the range of potential dynamics that this network exhibits under different cellular
T Gedeon: Preprint submitted to Elsevier Page 2 of 16
Dynamics of networks
conditions and different initial data; mathematicians can address these questions using differential equations only if
they are given parameters, nonlinearities and initial conditions. As discussed above, none of these are readily available.
Definition 1. A regulatory networkRN = (V ,E) is a graph with network nodes V = {1, 2,… , n} and signed, directed
edges E ⊂ V × V × {→, ⊣}. For i, j ∈ V , we will use the notation (i, j) ∈ E to denote a directed edge from i to j
of either sign; i → j will denote an activation or positive interaction, and i ⊣ j will denote a repression or negative
interaction.
We define the targets of a node i as
T(i) ∶= {j ∣ (i, j) ∈ E}
and the sources of a node i as
S(i) ∶= {j ∣ (j, i) ∈ E}
The goal of this paper is to describe two approaches that find a way to move beyond traditional ODE models
for network dynamics. Several ideas converge together to combine the best features of the continuous and discrete
dynamics into a powerful new tool.
3. Boolean networks
A Boolean model of network’s dynamics starts with set B ∶= {0, 1}. A state of a network is an element x ∈ Bn.
A Boolean network is a map f = (f1, f2,… , fn) with f i ∶ Bn → B and
f ∶ Bn → Bn.
The dynamical system generated by iterates of f is called synchronous Boolean network. LetH be the n-dimensional
hypercube graph with vertices v ∈ Bn and edges between x, y ∈ Bn if, and only if, these states differ in exactly one
component. Each synchronous Boolean map f induces asynchronous dynamics on H , given by a multi-valued map
F ∶ H ⇉ H , encoded by a state transition graph (STG), defined as follows: connect x ∈ H to all its neighbors y in
H with |xi − yi| = 1 for all indices i with the property that f (x) ≠ x .Boolean networks have been used to model gene regulatoiry netwoi rks (circuits) for many years [67, 63, 68, 57, 4,
18, 53]. Boolean network is determined by specification of a collection Boolean functions (f , f ,… , f ), which, in
turn, are each given by a truth table with number of inputs equal to number of incoming edge1s to2 the conrresponding
node in V . Any parameterization that encodes all collections of truth tables that are compatible with the network
RN will give parameterization of all Boolean models compatible with RN. However, since a change in a Boolean
network involves changing a value of some function f by a discrete unit (from 0 to 1, or vice versa) there is no natural
bifurcation theory that can delineate the region of valiidity of certain dynamics under changing cellular conditions.
In addition, synchronous dynamics also assume, unrealistically, that multiple genes change their states simulta-
neously. This problem is largely addressed by considering an asynchronous update rule for Boolean models. Lastly,
Boolean networks lack the means to ask questions about the time duration of other dynamic phenomena, like length
of the period of a cell cycle oscillation, since time is not a part of the modeling framework.
4. Continuous time models
In attempt to address some of these short-comings L.Glass [35, 36] introduced continuous time models with dis-
continuous right hand side, that are motivated by Boolean networks. They consider a system of differential equations
(see [28])
ẋi = −
xi + Fi(x̃1,… , x̃n) (1)
where
x̃i = 1 if xi > i; x̃i = 0 if xi < i. (2)
The system (1) is sometimes referred to as a Glass system. The Glass models have been used extensively [69, 45, 46,
6, 9, 50, 10, 15] in models of molecular network dynamics.
T Gedeon: Preprint submitted to Elsevier Page 3 of 16
Dynamics of networks
The connection to Boolean networks is straightforward; The function F = (f1,… , fn) reduces to a Booleanfunction f above, if we restrict the range of F to Bn. Note that the values of F ∈ Rn do not depend on values of x
but only on the values of x̃; therefore the range of F is finite - it can have at most 2n values. More formally, Rn is
decomposed into 2n disjoint open sets, parameterized by  ∈ Bn,
D() ∶= {x ∈ Rn | xi < i when 1 = 0, and xi > i when 1 = 1}.
The solutions of (1) inD() converge toward F () , but only until they leave the domain D(). Connecting trajectories
in neighboring domains that share a threshold, w
 henever possible, results in a piece-wise linear flow. There are many
subtleties as to when such merging of solutions in neighboring domains is possible, as well as how to extend these
solutions if it is not possible. We address these issues briefly below.
In Glass networks the function fi, viewed as a function Bn → B could be an arbitrary Boolean function. A moreappropriate class of Boolean functions for modeling biological networks is the class of monotone (unate) Boolean
functions [1, 21]
Definition 2. A function f ∶ Bn → B is increasing (decreasing) in x̃i, i ∈ {1,… , n}, f|x̃ =0 ≤ f ( f ≤ fi |x̃i=1 |xi=1 |x̃i=0
). In either case, f is said to be monotone (or unate ) in x̃i. It is said to be monotone, if it is monotone in all x̃i. In
particular, if it is increasing in all x̃i then it is positive.
For example, both logical functions AND and OR are monotone (increasing). However, the function XOR is not mono-
tone. The classes of monotone and positive Boolean functions attracted the interests of different research communities.
For review of their properties, we refer reader to Crama and Hammer [21] and Korshunov [1].
The concept of amonotone function captures the information embedded in the gene regulatory network where each
edge is either marked as activating or repressing [27, 70]. The number of monotone Boolean functions on n variables
is much smaller than that of the general Boolean functions, but it is still very large. In fact, the precise number of
monotone Boolean functions is the Dedekind number which is not known for n ≥ 9.
A different formalism that incorporates monotonicity implied by signed interactions in a regulatory network RN
was used by R. Thomas and coworkers [65, 66, 68, 67, 60, 61, 64]. They proposed that if a node i affects k downstream
nodes, then generically there should be k different thresholds ji at which the effect of i on a downstream node j isactivated.
The form of the differential equation (1) remains the same, but the definition of functionF changes. In Snoussi [60],
Fi consists of positive combinations of products and sums of contributions
i
±(x ) (see definition below) from nodes
j ∈ S(i). In a more general framework [61, 64, 56] the contributions are summedj
∑
Fi = ±(i, !), (3)
!⊂S(i)
but the sum ranges, in general, over all subsets of the input set of nodes S(i) of i [63, 56]. The form of functions ±(xj)is
{ {
+(i, !) = ki,! if xj > ij , − ki,! if xj < ij ,ij 0 if xj < ij , ij(i, !) = 0 if xj >  .ij ,
For every node i ∈ V the set of source nodes of i decomposes S(i) ∶= S+(i) ∪ S−(i) into the set of positive and
negative regulators of i respectively, by requiring that j ∈ S+(j) (j ∈ S−(j)) if and only if j → i (j ⊣ i). The
collection of numbers K ∶= {ki,!} satisfies the monotonicity assumption
k ′i,! ≤ ki,!′ for ! ⊂ !
for subsets ! ⊂ !′ ⊂ S+(i) of positive inputs. The monotonicity condition reverses
ki,! ≥ ki,!′ for ! ⊂ !′
for for subsets ! ⊂ !′ ⊂ S−(i), see [11].
The motivating assumption is that the basic expression level of j-th input is k for ! = {j}; the simultaneous
expression of inputs j and s, that corresponds to ! = {j, s}, will be at the level i,jki,js that is larger than both ki,j and
ki,s. We call the system (1) with nonlinearities (3) a monotone system.
T Gedeon: Preprint submitted to Elsevier Page 4 of 16
Dynamics of networks
There are other slightly different formulations of this system. For instance Ironi [43, 44] replaces x̃j in generalformulation (1) by a Boolean function with values in B, i.e. Z (x ,  ) ∈ B, and assumes that for any j the thresholds
ij of variable xj are distinct and that F are multilinear polyn
iojmijalsijin the Z .
Chronologically the most recent formi ulation of the ODE system with disciojntinuous right hand side, which gener-
alizes Glass system (1), was considered in [26]. We will call this a switching system. This formulation was used in a
comprehensive description of the dynamics across parameter space, called DSGRN (Dynamic Signatures Generated
by Regulatory Networks), which is presented in next section.
Consider
ẋi = −
ixi + Λi(
±
i1(x1),… , 
±
in(xn)) (4)
where each entry is either function + or a function −, defined by
{ {
+ij(xj) =
Uij if xj > ij , −
L if x <  , ij(x ) =
Uij if xif j
< ij ,
j
ij j ij Lij x
.
j > ij ,
Assume that the functions Λi(y1,… , yn), i = 1,… , n are of the form
∏∑
(yj + yj +…+ yji1 i2 ik(j)),j
where each ys occurs exactly once. Such a function is automatically multi-linear i.e. linear with respect to each y ,and all coefficients are 1. The values Lij < Uij represent lower and upper activation (repression) level of gene
i
i by
gene j. Note that +ij mediates an up-regulation of i by j, while − mediates a down-regulation.We describe the difference between monotone system (3) andijswitching system (4). The first difference is that the
nonlinearity Fi is a sum in (3) but function Λi in (4) is arbitrary multi-linear function. This provides more freedom tothe switching system to match the choice of nonlinearity to biology represented by the networkRN.On the other hand,
the collection of constants K for (3) is only required to meet monotonicity assumption, but is otherwise arbitrary. The
collectionK can be computed for any choice of multilinear function Λi, but the type of multilinear function that has tobe selected for each i in switching system (4) will restrict in subtle way what collectionsK are attainable as parameters
Lij < Uij are allowed to vary.We now describe the differences between Glass system (1) and switching system (4). First, (4) allows decay
constants i to vary between nodes i. More importantly, there is a difference between how the parameters enter function
Fi vs. Λi. In (1) the coefficients of function Fi parameterize expression level of gene x , while in (4) the expressionlevel of xi is determined directly by the values of
i
U ,L
how they are combined to the final value of . Rijestriij
along the input edges to i, together with Λi that determines
Λ cting the range of each function ± to U = 1 and L = 0,
the function + ( or −) becomes an increasingi (decreasing) Boolean function of the Boolean variable x̃ (see (2)).
Combining such functions inΛ results in a monotone Boolean function. Using general valuesL,U instead oj f Boolean
values 0, 1 represents one way iin which the switching system generalizes the Glass model.
The second difference comes from the number of different thresholds that (4) considers. In Glass system (1) when
xi > i then the state x̃i of gene i is x̃i = 1 and this value enters the function F (x̃). In other words, all the downstreamnodes are activated. On the other hand, in a switching system (4) a node i activates downstream nodes xj at their ownthresholds ji and it may activate some, but not others of the downstream nodes. This depends on. the relative size ofthe value of the activated gene Λi and the thresholds associated to each downstream edge under a generic condition
ji ≠ ki for all i and all k ≠ j. (5)
While the continuous time systems (1,3,4) were developed to address some deficiencies of Boolean networks related
to the lack of continuity of time, space and parameters, these systems generate their own sets of technical problems.
The main technical difficulty in the traditional view of (1,3, 4) when viewed as systems of ordinary differential
equations, is in merging solutions that exist in each domain (Figure 1a) to global solutions that are defined for all t > 0.
One problem is posed by the solutions in domainD that enter co-dimension k-manifolds with k ≥ 2 that are formed as
intersections of more than two neighboring domains. The more serious problem that affects an open set of solutions
is when the solutions from neighboring domains both enter their common co-dimension one intersection which is a
T Gedeon: Preprint submitted to Elsevier Page 5 of 16
Dynamics of networks
ji xi ii xi
(a) Solutions can be joined across a threshold ji atwhich x affects x . (b) Negative self-regulation of node i at threshold ii.i j
Figure 1:
subset of {x ∈ ℝn | xi = ji} (Figure 1b.) It is easy to see that this happens if, and only if, the node i ∈ RN has anegative self-loop.
There has been a substantial body of work devoted to the proper definition of solutions of (1,3, 4) for all t >
0. They either invoke Filippov’s idea of differential inclusions and replacing (1,3, 4) by corresponding differential
inclusions [31], or consider only those solutions that are limits of solutions of nearby steep sigmoidal functions ([43,
44, 71]. However, these issues distract from the fundamental problem of developing a modeling toolbox for biology
and often form a un-surmountable block for biologists who may otherwise be tempted to use it.
It is our view that the system (4) should be viewed as a computational model for computing lattices of attracting
blocks, which are sets in the phase space that are positively invariant and form neighborhoods of the attractors. In can
be shown, that the lattice of attracting blocks of (4) is part of the lattice of attracting blocks for all neighboring smooth
ODE systems [34]. This shifts focus from the individual trajectories to the attracting sets, which are robust under the
change of nonlinearity from piecewise constant to smooth.
In order to avoid this technical difficulty, the following assumption is often made
Assumption: The network RN has no negative self-regulation.
In the context of the applications, a network is usually a simplification of a more complex network and inclusion
of an intermediary like an mRNA, or another protein into a negative self-loop allows this assumption to be satisfied.
4.1. State transition graph
Both Boolean maps and ordinary differential equation (ODE) systems (1, 3,4) admit a discrete description of the
phase space. For Boolean map f , we will always represent dynamics by iterates of its multivalued asynchronous
update function F ∶ H ⇉ H , which can be viewed as a graph on n-dimensional hypercubeH .
For the Glass ODE system (1) the underlying phase space for differential equations is ℝn. However, the phase
space Rn decomposes to 2n domains D() where  ∈ Bn and there are bijections between the set of domains D(),
the set  ∈ Bn, and the set of vertices ofH .
The state transition graph (STG) for ODE models (1) on H is defined analogously to the Boolean model. We
connect  ∈ H to all its neighbors ′ inH with ′|i − i | = 1 for all indices i with the property that
F ()
[ ]i() ≠ 
 i
.
The STG can be viewed as a multi-valued map FG ∶ H ⇉ H , analogous to an asynchronous update of the Booleanfunction F ∶ H ⇉ H .
The connection between the dynamics of iterates of F ∶ H ⇉ H and the solutions of ODE (1) is the at the heart
of the approach described in this paper. Its origin is the oGbservation that all ODE solutions in D() converge toward
F ()

 , while in the domain D. The discrete description of the continuous dynamics has a natural appeal, but the detailsof how exactly the long term behavior of the solutions of ODE correspond to the long term behavior of iterates of F
is subtle [43, 44, 50, 71, 28]. Any attempt to make this comparison on the level of an individual solution must addresGs
the technical issues outlined above.
T Gedeon: Preprint submitted to Elsevier Page 6 of 16
Dynamics of networks
For more general ODE systems (3, 4) the key assumption (5) implies that each variable x has |T (i)| thresholds
(recall that T (i) is a set of out edges from a node i in the network RN). Therefore, by analogy to isystem (1) where each
variable has one threshold, the continuous phase space Rn consists of Πni=1(|T (i)| + 1) domains D() where  is nowa multi-index with i ∈ {0, 1, 2,… , |T (i)|}. Let
H̄ = { |  = (1,… , n), i ∈ {0, 1, 2,… , |T (i)|}}
be a generalized hypercube consisting of nodes with labels . Each node represents a corresponding domain D() ⊂
Rn.
We construct a state transition graph (STG) on H̄ using the fact that solutions in each domainD() converge toward
its associated target point T ()) = (Λ ()∕
 ,… , (Λ ()∕
 ). STG is a graph representation of asynchronous update
of the discrete-valued function 1 1 n n
Λ ∶ D()→ T ().
This means that STG is a graph on H̄ with edges between nodes ,  with | − | ≤ 1 assigned in the following way
• if T () ∈ D() then there is a self-edge  → ;
• if | − | = 1, i = i − 1 and the threshold between D(), D() is , then
1.  →  if Ti() > ;
2.  →  if Ti() < .
Note that the the latter assignment is well defined unless both conditions are satisfied at the same time. Then it must be
that the threshold  is a threshold of variable x that affects the function Λ , which implies that node i regulates itself
n the network RN. Furthermore, this self reguliation must be negative, sincie the Λ value decreases when variable x
increases from domainD() to domainD(). Since we assume that there are no negi ative self-loops in RN, the STG isi
well defined and there is at most one edge between neighboring nodes in H̄ . STG can be also viewed as multi-valued
map  ∶ H̄ → H̄ .
The state transition graph is a finite object that can be interrogated in different ways. The Boolean community has
implemented a variety of algorithms that use a correspondence between the states of the STG and their description by
logical expressions. This resulted in very fast algorithms that decide reachability of certain states from other states;
we only mention the selected work on software implementation [30, 38] and applications [12, 5].
The DSGRN approach [26, 25, 24, 33] focuses on the compression of the information about the dynamics of  by
computing strongly connected path components of STG and the reachability properties between them. The strongly
connected path components are nodes of a Morse graph which is a Haase diagram of a partial order imposed by the
reachability on the nodes of Morse graph, that is inherited from reachability in STG. Morse graph is an analogous
concept to that of a Morse decomposition [20], restricted to a finite dynamical system  ∶ H̄ → H̄ .
Each attractor A of  i.e. a positively invariant collection of nodes in H̄ under  , corresponds to a set of domains
D() ⊂ Rn that is positively invariant under the flow of the switching system. This set of domainsmay consist of several
disjoint non-contiguous collections of domains. The union of these domains is called an attracting block corresponding
to A. The boundary of such an attracting block is transverse to the solutions of (4) and are therefore robust under and
C1 perturbations of the flow, and therefore C0 perturbation to the nonlinearities. This class of perturbations includes
smooth systems in the neighborhood of the piece-wise constant nonlinearities. These ideas have been used in [34] to
show that in a two dimensional system under mild conditions, there is an injective map respecting the partial order that
maps the Morse graph of the switching system (4) to a Morse decomposition of any differential smooth equation with
nonlinearities C0-close to switching system (4). This shows that the switching system provides a rigorous lower bound
on complexity of dynamics of nearby smooth systems. Since switching systems cannot directly capture the existence
of singular solutions [68, 63, 56], the fact that switching system only provides a lower bound is not surprising.
5. Dynamics in context of changing parameters.
We now come to the main part of the paper where we discuss how to study either Boolean maps or ODE systems
(1-4) as the Boolean functions or nonlinearities change. In the ODE systems, this is called bifurcation theory where
T Gedeon: Preprint submitted to Elsevier Page 7 of 16
Dynamics of networks
one seeks to find parameter values where qualitative features of the dynamics abruptly change. Bifurcation values of
the parameter bound areas of the parameter space where the dynamics is qualitatively the same.
It has been recognized at least as early as [63] that by systematic change in parameters k in system (3) one obtains
different qualitative dynamics in STG. This has been formalized in [56, 11] by introducinig,!the concept of qualitative
parameters, which parameterize potentially different dynamics of the state transition graphs. The recent work [2],
inspired by [63], uses the correspondence between Boolean maps and Glass systems (1) to define bifurcation theory
for Boolean maps. Recall that for a Boolean function f ∶ Bn → B one can define the set of truth values T (f ) and
false values F (f ) by
T (f ) ∶= { ∈ Bn | f () = 1}, F (f ) ∶= { ∈ Bn | f () = 0}.
Using this concept, a very recent work [27] organizes the structure of the set of monotone Boolean functions f by
ordering the sets T (f ) by inclusion. i
A natural parameiterization for systems (4) has 3E + N continuous parameters (L ,U , and  for each edge
(j, i) and 
i for each node i). We will describe the relationship between this continuous p
ijaramijeter spaicje and the finite
number of admissible state transition graphs. We then compare this description to the other approaches.
5.1. Parameter graph
The parameter space of system (4) consists of a collection
Q ∶= {(Lij , Uij , ij , 
i) | for all (i, j) ∈ E, i ∈ V , such that Lij ≤ Uij .} (6)
We will consider an open subset P ⊂ Q where
Λi ≠ 
iij for all i, j.
We make a following definition that leads to a concise description of all dynamics that can be generated by a given
network across all parameters.
Definition 3. For node i let Zi be the set of input strings i = (a1,… , ap ) of the length pi ∶= |S(i)|, the number ofi
sources of i, where each aj ∈ {Lij , Uij}. A state of a node i at a parameter p ∈ P is a collection of inequalities
Σi(p) ∶= {Λi(i) < 
iji or Λi(i) > 
iji | j ∈ T(i), i ∈ Zi}
where for each pair (i, j) there is exactly one inequality (< or>) in the collection Σi. Importantly, sinceΛi is assumed
to be a product of sums where each dependent variable occurs at most once, Λi is monotone in each of its arguments.
The choice of inequalities has to respect this monotonicity.
A state of the system is a collection Σ(p) = (Σ1(p),Σ2(p),… ,Σn(p)).
We make three key observations [34]
1. The parameter p ∈ P , and the identity of domainD() determines uniquely the values of strings  ∶= (1(),… , n()).
2. Given  , the state of the systemΣ(p) determines position of target points of all domainsD() and thus determines
unique state transition graph.
3. The number of states is finite, and so there are open sets of parameters p ∈ P that generate the same state Σ(p).
A parameter graph PG for a network RN has a node for each state of the system Σ; an undirected edge between
two states Σ and Σ exists, if Σ and Σ differ by exactly one inequality.
Because1the stat2es of the syst1em dep2end on the algebraic form of Λ , the parameter graph will be different for a dif-
ferent collection of algebraic forms (Λ ,… ,Λ ). Each assumed algebriaic form of the function Λ imposes constraints
on which the monotone Boolean funct1ions arenrealizable by that algebraic form [26]. i
It is worth noting that the set of parameter nodes is equivalent to the set of qualitative parameters discussed in
Richard et. al. [56] and Bernot et. al. [11]. The parameter graph PG organizes this set into a graph that captures
geometrical structure of the corresponding parameter regions in the parameter space.
Note that since each state Σ can be written as Σ = (Σ1,… ,Σn), the parameter graph PG is a product of parametergraphs for each node i ∈ RN
PG = Πni=1PG(i).
T Gedeon: Preprint submitted to Elsevier Page 8 of 16
Dynamics of networks
The product structure significantly simplifies construction of PG, since PG(i) only depends on the type of the node
i; i.e the number of input edges, the number of output edges and the algebraic form of function Λ that combines the
inputs. These graphs can be precomputed and then combined together to construct PG for any netwi ork RN.
In Figure 2(a) we show PG for a network node z with two inputs from nodes x, y and one output. Here we assume
that Λ(+(x), +(y)) where x, y are the continuous variables representing the concentration of chemical species and
Λ(u, v) = uv is a product. In order to show connection between DSGRN parameter graph and qualitative parameters
of [56, 11], we associate logical variables X to x and Y to y
{ {
X ∶= 1 x > zx0 x <  Y ∶=
1 x > zy
zx 0 x < zy
which represent Booleanization of +(x), +(y), respectively. In Figure 2(a) we use subscriptsX and Y inL ,U ,L , U
instead of Lzx, Uzx, Lyx, Uyx to enable comparison to Boolean description below.
X X Y Y
In order to present a succinct representation of parameter nodes, we only include in Figure 2(a) the minimal set of
inequalities from which the other inequalities characterizing the parameter node that can be deduced by using standard
inequalitiesLX < UX , LY < UY . For instance, the top node in Figure 2(a) also includes implicit inequalitiesL L <
UXLY , LXLY < LXUY , UXLY < UXUY , LXUY < UXUY . Observe that each node differs from neighboring
XnoYdes
by exactly one inequality. We also note that the structure of this parameter graph would be the same if we chose
Λ(u, v) = u + v. For functions Λ with more that two inputs, the parameter graphs may be different for different
algebraic forms of Λ.
We now connect the parameter graph to the language of Booleanmaps. Consider a parameter nodeΣ = (Σ ,… ,Σ )
and assume that the j node of RN has k-inputs and one output node 1|S(j)| = k, |T(j)| = 1. We demonstrate thnat
the value Σj corresponds to unique Boolean function g ∶ Bk → B. As a first step, we map every input sequence
j = (a1,… , ak), ai ∈ {Lji, Uji} to a sequence  ∈ Bkj by mapping each U to 1 and L to 0. If an input sequence
j ∈ {L,U}k is such that the inequality Λj(j) < 
 belongs to Σj , then we assign the corresponding sequence
 kj ∈ B to the false set j ∈ F (g). Similarly, if an input sequence j is such that the inequality Λj(j) > 
 belongs
to Σj , then we assign the corresponding sequence  ∈ Bkj to the true set j ∈ T (g). This uniquely associates to
Σj a Boolean function g for which F = F (g) and T = T (g) with Bk = T (g) ∪ F (g). The collection {T (g)}g over
all Boolean functions g ∈ Bk → B, ordered by inclusion, produces a partial order that recovers the lattice structure
of Figure 2(a). Each set of true values T (g) gives rise to a description of the function g in terms of logical formulas
depicted in Figure 2(b). Therefore, there is one-to-one correspondence between parameter graph in Figure 2(a) and
the graph of Boolean functions Figure 2(b).
Consider now a parameter node Σ = (Σ1,… ,Σn)where the j node ofRN has k-inputs and l output nodes |S(j)| =
k, |T(j)| = l with corresponding thresholds 1 <  < … <  . The value Σ corresponds to a collection of lmonotone Boolean functions g ∶= (g1, g2,… , gl) wh
2ere the set olf true values jT (gi) of a function g is a collectionof  ∈ Bk for which the corresponding  satisfy Λ ( ) > 
  . The ordering of the thresholds impliies that the true
valujes of the collection g satisfy j j j j i
T (gl) ⊆ T (gl−1) ⊆… T (g1).
In the Figure 2(c) is the PG(i) for a network node i with one input and two output edges. Since each output edge
corresponds to a threshold, there are two thresholds 1 and 2. To explain the shape of this graph, assume first that
1 < 2. If L < U denote two possible values of the input, then we have the following collection of parameter nodes
(a1) U < 1, (b1) L < 1 < U < 2, (c1) 1 < L < U > 2;
(d1) L < 1 < 2 < U ; (e1) 1 < L < 2 < U ; (f1) 2 < L
The nodes a2,… , f2 in Figure 2(c) with the subscript 2 are obtained by exchanging  and  in the above formulas,which is the equivalent to changing their mutual order to 2 < 1. Nodes that differ by a
1single i2nequality are connected
by an edge.
Computations of parameter graph PG are limited by memory storage, rather than the computational complexity.
DSGRN databases with billion (109) nodes has been computed. These typically correspond to networks with 6-10
nodes. The precise number depends more sensitively on the complexity of regulation at each node, rather than on the
T Gedeon: Preprint submitted to Elsevier Page 9 of 16
Dynamics of networks

 < LXLY
LXLY < 
 < {UXLY , LXUY } True f1 f2
e e
UXLY < 
 < LXUY LXUY < 
 < U L
X ∨ Y 1 2
X Y
Y X c1 d1 d2 c2
{UXLY , LXUY } < 
 < UXUY
X ∧ Y b1 b2
UXUY < 

False a1 a2
(a) (b) (c)
Figure 2: (a) Parameter graph for two input edges from nodes x and y with inputs LX < UX and LY < UY respectively,
and one output edge. Each node is marked by the set of inequalities; we only include minimal set of inequalities sufficient
to reconstruct all of them. Nodes marked red (thicker boxes) are essential parameter nodes. (b) Description of the
corresponding monotone Boolean function; (c) Parameter graph for one input and two output edges. Symmetry between
left and right corresponds to change of order of output thresholds.
number of nodes themselves. If a node has only up to two inputs and outputs, the parameter graphs is rather small
O(101); number of inputs increases this number extremely rapidly.
One of the way to reduce the size of the parameter graph is to employ the concept of an essential parameter
node [26, 24, 22], related to the concept of an essential monotone Boolean function [27, 70]. A parameter node
Σ = (Σ1,… ,Σn) is essential if for every Σj , j = 1,… , n the following holds:
a For each threshold ij associated to i ∈ T(j), there is a input string  such that Λj( ) < 
jij and an input string ̄
such that Λj(̄ ) > 
jij ;
b For every k ∈ S(j) there exist states 1, 2 that differ only in k−th component, and a threshold ij such that
Λ(1) < 
jij < Λ(2) or Λ(2) < 
jij < Λ(1).
In other words, each input edge and output edge plays a role in determining an output of the network node j.
The concept of the essential monotone function [27, 70] coincides with the concept of the essential parameter node
for network nodes with a single output. In Figure 2(b) the essential Boolean functions are in red and correspond to
those functions where both inputs X and Y play a role in determining the outcome.
A parameter node is non-essential if it is not essential. Non-essential parameter nodes can be viewed as param-
eter nodes of a reduced network where at least one of the input, or output edges, can be removed from the network
(Figure 2(b)). DSGRN database restricted to essential parameter nodes is smaller and thus potentially computable for
larger networks.
Definition 4. A DSGRN (Database of Signatures Generated by Regulatory Networks) for a network RN
⋃
MG(p),
p∈PG
is a collection of Morse graphs parameterized by the nodes of the parameter graph PG.
The idea of considering a class of nonlinearities Λ in (4) (or F in (1)), rather than analyzing a single Boolean or
switching network was considered by several authors. Iit has been reicognized by [29] that precise parameter values can
T Gedeon: Preprint submitted to Elsevier Page 10 of 16
Dynamics of networks
be replaced by inequalities between parameters to get an identical qualitative description of the dynamics in terms of
state transitions. This has been used in [43, 7] to iteratively generate inequalities between the parameters at which the
system (4) qualitatively matches the time series of the experimental data by state transitions in STG. In the language
of DSGRN, as more of the time series is known, a smaller number of parameter nodes can reproduce the experimental
time series. This iterative process zooms in on the parameter regions of interest. These papers use greedy methods
which may not find all regions of interest. Compared to DSGRN, the advantage of these methods is, that since it does
not compute the entire parameter space, it may be able to handle larger networks; the disadvantage is that for every
new time series, the computation needs to be repeated from the start.
More recently [2] introduced logical bifurcation diagrams . This is a way to track changes in the state transition
graph as a function of a particular parameter. In the language of this paper, this corresponds to a restriction of DSGRN
to a particular subgraph PG(i) for some network node i and considering changes in Morse graph as a function of paths
through parameter nodes. The paths of interest are either increasing or decreasing paths in the graphs PG(i) like those
in Figure 2(a) (c). Since the entry point of [2] is the Boolean network formalism, which is discrete, they introduce
the ODE system (1) to embed Boolean functions into a setting with varying parameters in the way that is described
above. The emphasis is on the generation of one-dimensional paths through parameter graphs starting at a particular
parameter node in PG(i) Since PG(i) is a partially ordered set there can be multiple choices for these paths.
The work [27] shows how to traverse a graph that only includes essential monotone Boolean functions with k
inputs. This can generate potential bifurcation paths within essential monotone Boolean functions and can be used to
define a neighborhood of an essential monotone Boolean function within this class.
6. Applications
The advantage of the description of the dynamics by DSGRN database is the ability to ask questions that cannot be
asked by considering a single Boolean function or an ODE at a single parameter value. Global information about the
dynamics of the system is represented in DSGRN in different ways. For instance, the nodes of the Morse graph that
have no outgoing edges represent invariant sets that are stable in the sense that there are no other Morse nodes that can
be reached from this node. We will call such nodes stable. If a Morse graph has two such minimal nodes, the system
exhibits bistability; more nodes indicate multi-stability.
Furthermore, we use annotation of the Morse nodes of Morse graph at p ∈ PG as a way to represent a dynamical
signature of the dynamics of all parameters that belong to the parameter region represented by p. Currently we use
three types of annotation. A Morse node is designated FP if it corresponds to a single node  ∈ H̄ and there is a
self-loop  →  in STG. The name FP (for "fixed point") signifies that the corresponding domain D() ⊂ Rn will
contain a fixed point of the dynamics of (1,3) or (4). Annotation XC (for Partial Cycle) designates a Morse node
that represents a set of domains D(i) where there is at least one component j that is constant across the collection
i, i.e. [i]j = const for all i. In other words, such a component does not change its state in the set of domains
D(i). Annotation FC (for Full Cycle) describes a Morse node that represents a set of domains D( ) where noneof the components are constant across all i. We hasten to add that the existence of a stable Morse nodie FC does not
guarantee that for every parameter corresponding to this parameter node, the ODE (4) has a stable periodic orbit, as has
been shown by Glass [37]. However, we view such a set of periodic domains as an adequate proxy for the periodicity
observed in the experimental data, which is both finite and subject to experimental noise.
One way to use annotated Morse graphs in DSGRN database is to rank networks based on their ability to robustly
exhibit an oscillatory dynamics, based on prevalence of stable FC in the DSGRN database. This may be used to rank
potential networks that may explain oscillatory phenomena like a cell cycle, a circadian rhythm, or the emergence of
a malaria parasite from red blood cells. A network that does not show any FC in the entire parameter graph is clearly
not an appropriate model for a system that exhibits periodic behavior. This trivial conclusion may lead to the rejection
of many potential networks, and thus may lead to a substantial reduction of hypothesis space. This rejection can be
enhanced by assuming that the network must robustly exhibit periodicity represented by FC; this allows a ranking of
the networks based on the FC prevalence and a rejection of lower ranked networks.
An additional restriction comes from considering extrema (i.e minima and maxima) of an experimental time series.
A path in a state transition graph of a switching system (4) also corresponds to an ordered set of minima and maxima
of the coordinate functions in (4). Therefore we can evaluate not only how prevalent is the periodic behavior encoded
in Morse set FC, but by comparing the sequence of maxima and minima along an arbitrary path through such an FC
with the ordered set of minima and maxima of the same genes in experimental data, we can reject those networks that
T Gedeon: Preprint submitted to Elsevier Page 11 of 16
Dynamics of networks
do not match the experimental ordering or extrema [25],
We describe a broad set of questions that can be addressed by the approaches described in this article.
6.1. Delineation of network dynamics
E2F-Rb network [48, 3, 54, 16] plays a key role in one of the most important decisions in the human cell life
cycle: whether to start DNA replication and initiate proliferation. This decision is based on several factors, but once
the process has started DNA replication, it must finish. Therefore, the influence of the factors that led to the decision
must be uncoupled at the moment of the decision, called the restriction point of the cell cycle [52, 14, 59]. Not
surprisingly, mutations in E2F-Rb network are linked to many cancers [17, 19, 49]. The requirement of irreversibility
and decoupling suggests that phenotypically a bistable switch may underlie the restriction point.
In [33] the authors analyzed a series of networks inspired by the work of Yao and his collaborators [76, 77]. In
([77]) the authors simplified a large E2F-Rb network by clustering nodes into three groups, with a number of potential
edges between them that were consistent with the edges of the larger network. Then they numerically simulated many
potential 3-node networks consistent with this information at a large, but fixed, set of parameter values, and evaluated
which of these networks most robustly exhibit different types of hysteresis. In [33] this analysis was first repeated by
computingDSGRN for all 3-node compatible networks. Then the authors went a step further and looked at eight 6-node
networks that were expanded from the 3-node networks and asked which of them most robustly showed the hysteresis
behavior. Interestingly, the network that exhibits hysteresis most robustly according to the DSGRN database, matched
the topology of the network that performed the same function in yeast. Since the genes involved in the restriction point
in yeast share no homology with the genes in humans [23], this result suggested that the evolution may have selected
the same network topology independently multiple times to perform the same function.
In order to search for hysteresis in the E2F-Rb, [33] searched paths in the parameter graph PG(S) where S was
a node in the graph affected by the growth rate of the cell. Along every path they asked if the set of Morse graphs
exhibited the following sequence of events: first, a single FP that represents non-proliferative state, followed by aMorse
graph exhibiting bistability, and finally, at the end of the path, a Morse set with single FP representing a proliferative
state. This annotated Morse graph representation of hysteresis allowed [33] to enumerate prevalence of hysteresis in
parameter graphs of different networks.
6.2. Control and manipulation of dynamics
There have been several attempts to use predictions of Boolean networks to control and manipulate their dynamics.
The program kali developed by ([55]) considered a single Boolean network and assumed that the attractors, and in
particular the fixed point attractors, corresponded to network phenotypes. The goal of the program was to reduce the
size of the basins of attraction of the so called pathological attractors that correspond to un-desired outcome of the
network dynamics. This program allows a choice between the synchronous and asynchronous updates and allows for
multi-valued logic.
Still within a single Boolean network, the series of papers [78, 80, 79, 81] assumed that there was a set of con-
trol nodes that can be externally manipulated. They used optimal control techniques to derive a control policy that
minimized the likelihood of transitioning into an undesirable state in a stochastic context.
A slightly different approach was taken by Biane and Delaplace [13] who extended the Boolean formalism to
Boolean control network, where additional control Boolean nodes were added to the system. These nodes are not dy-
namically updated, but they serve as preset control variables. Including them in a logical expression with dynamically
updated variables of the Boolean network allows encoding of the node control into extended Boolean framework. One
can view this approach as a discrete parameterization of a family of Boolean maps by all possible Boolean values of
the control nodes. It is possible to translate the Boolean control framework into the framework of DSGRN.
A mathematically different approach was taken by ([51]). They translated the problem of finding control point
candidates to a problem of solving a set of polynomial equations. These were then solved by a powerful techniques of
computational algebra. The output was a set of actions that are capable to steer the network towards desirable states.
6.3. Learning parameters that fit the data
An important problem is to match a Boolean or an ODE model to the experimental data, which is often presented
in the form of a time series of several (or all) network variables.
As mentioned above, the work of ([7]) and ([43]) consider ODE systems (1)-(4) with the goal of selecting a set of
parameter nodes which produce solutions that match experimental data. The key observation is that the information
T Gedeon: Preprint submitted to Elsevier Page 12 of 16
Dynamics of networks
that needs to be matched is qualitative, since (1)-(4) produce piecewise linear solutions. Both papers ([7])-([43]) start
with the assumption that the experimental trajectory has been mapped into a sequence of domains D(). Then, given
an initial domain, ([43]) construct a tree where each branch connects to all neighboring domains and each edge is
annotated by a list of parameter inequalities that are consistent with that transition. This results in a very large graph
and is not practical for networks with more than a few nodes, even for a relatively short time series.
The paper ([7]) encodes qualitative observations about data into a temporal logic formula and uses model checkers
([5, 8, 12, 32]) to check if a state transition graph generated at a particular parameter node satisfies this formula. To
avoid the explicit construction of the parameter graph PG they use an alternative approach. They effectively use the
branching graph structure of ([43]) without explicitly constructing it. It this way, for a given path through the state
transition graph that is generated by the experimental data, the work ([7]) computes parameter constraints that are
consistent with the path. These constraints become more restrictive as the length of the path increases.
Another approach to the matching of an experimental time series to the dynamics of network has been done in the
context of DSGRN [25]. Given an experimental time series and a proposed network, the goal is to find all parameter
nodes of the parameter graph where the dynamics is consistent with the time series. To do this, a partial order of the
maxima and the minima of individual time series is extracted from the time series and compared to STG dynamics
generated by all parameter nodes of the DSGRN database. The consideration of a partial order rather than a full
order is motivated by the fact that experimental time series are often sparsely sampled and multiple genes can achieve
maximum or minimum at the same time point. The result of the analysis is a collection of parameter nodes where the
discrete STG dynamics exhibits the same partial order of maxima and minima as the experimental time series. Then
networks can be compared based on robustness of matching i.e. prevalence of this behavior across all parameters.
There has been a lot of work devoted to learning Boolean networks from data. In the context of a parameter graph
PG, this is equivalent to selecting which parameter nodes are consistent with the observed dynamics. The paper [39]
introduces caspo, a free open-source tool to learn the (Boolean) logic models of signal transduction. It uses CellNOpt
pre- and post-processing routines [62]. The output is a set of Boolean models that are capable of reproducing data. In
our language, the output is a collection of parameter nodes at which the state transition graph STG has a path that is
consistent with the data. The follow-up paper ([73]) goes a step further. Given a collection of potential Boolean models
that are consistent with the data, it proposes a minimal set of experiments that are capable of distinguishing between
them. In this case, however, the method is restricted only to Boolean networks representing RN with no feedback
loops.
Finally, ([72]) used methods of computational algebra to learn all Boolean maps that are consistent with the data.
The algorithm consists of encoding possible wiring diagrams using ideals and algebraic sets and choosing those that
are minimal using primary decomposition and irreducible components.
7. Discussion
Wedescribe how relatively independent ideas developed over the last 50 years in dynamical systems (Glass systems,
switching systems) and computer science (Boolean networks, automatic model checking), started to merge into a
unified modeling platform for systems biology. Facing the challenge of describing the dynamics of networks, where
the nonlinearities, parameters and initial data are known only qualitatively, a combination of continuous time dynamics
and its discrete time description allows enumeration of all potential dynamics compatible with the network structure.
This description, by necessity, matches a qualitative input that systems biology is able to provide. However, this
approach allows
• the delineation and classification of all dynamics consistent with network structure;
• the description of changes in dynamics induced by the changes in parameters;
• the evaluation of network structures for their potential fit to the experimental data;
• the rejection of proposed network structures that are incompatible with the experimental data.
These techniques, including perhaps the most clearly formulated approach called DSGRN [26, 25, 24, 34, 42, 22, 40],
have been successfully used in cell biology [33, 75].
We hope that this contribution will spur more applications, as well as further the development of these techniques.
T Gedeon: Preprint submitted to Elsevier Page 13 of 16
Dynamics of networks
References
[1] A, K., 2003. Monotone Boolean functions. Russ. Math. Surv. 58, 929–1001.
[2] Abou-Jaoude, W., Monteiro, P., 2019. On logical bifurcation diagrams. J. Theor. Biol. 466, 39–63.
[3] Adams, M., Sears, R., Nuckolls, F., Leone, G., Nevins, J., 2000. Complex transcriptional regulatory mechanisms control expression of the
e2f3 locus. Mol. Cell. Biol. 20, 3633Ð3639.
[4] Albert, R., Collins, J., Glass, L., 2013. Introduction to Focus Issue: Quantitative approaches to genetic networks. Chaos 23, 025001.
[5] Barnat, J., Briml, L., Cerna, I., Drazan, S., Fabrikova, J., Safranek, D., 2009. On algorithmic analysis of transcriptional regulation by ltl model
checking. Theor. Comput. Sci. 410, 3128–3148.
[6] Batt, G., Belta, C., Weiss, R., 2007a. Model checking genetic regulatory networks with parameter uncertainty, in: Bemporad, A., Bicchi, A.,
Buttazzo, G. (Eds.), Hybrid Systems: Computation and Control, Springer Berlin Heidelberg, Berlin, Heidelberg. pp. 61–75.
[7] Batt, G., Page, M., Cantone, I., Goessler, G., Monteiro, P., DeJong, H., 2010. Efficient parameter search for qualitative models of regulatory
networks using symbolic model checking. Bioinformatics 26, i603–i610.
[8] Batt, G., Ropers, D., de Jong, H., Geiselmann, J., Mateescu, R., Page, M., Schneider, D., 2005. Validation of qualitative models of genetic
regulatory networks by model checking: analysis of the nutritional stress response in escherichia coli. Bioinformatics 21, i19–28. doi:10.
1093/bioinformatics/bti1048.
[9] Batt, G., Yordanov, B., Weiss, R., Belta, C., 2007b. Robustness analysis and tuning of synthetic gene networks. Bioinformatics 23, 2415–2422.
[10] Bernard, O., Gouze, J., 2002. Global qualitative description of a class of nonlinear dynamical systems. Artificial Intelligence 136, 29–59.
doi:10.1016/S0004-3702(01)00169-2.
[11] Bernot, G., Cassez, F., Comet, J., Delaplace, F., Muller, C., Roux, O., 2004a. Application of formal methods to biological regulatory networks.
J. Theor. Biol 229, 339–348.
[12] Bernot, G., Comet, J., Richard, A., Guespin, J., 2004b. Application of formal methods to biological regulatory networks. J. Theor. Biol 229,
339–348.
[13] Biane, C., Delaplace, F., 2018. Causal reasoning on boolean control networks based on abduction: theory and application to cancer drug
discovery. IEEE/ACM Transactions on Computational Biology and Bioinformatics doi:10.1109/TCBB.2018.2889102.
[14] Blagosklonny, M., Pardee, A., 2002. The restriction point of the cell cycle. Cell Cycle 2, 102–109.
[15] Bornholt, S., 2008. Boolean network models of cellular regulation: prospects and limitations. J. R. Soc. Interface 2008 5, 134–150.
[16] Bouchard, C. qnd Thieke, K., Maier, A., Saffrich, R., Hanley-Hyde, J., Ansorge, W., Reed, S., Sicinski, P., Bartek, J., Eilers, M., 1999. Direct
induction of cyclin d2 by myc contributes to cell cycle progression and sequestration of p27. embo j. 18: 5321-5333. EMBO J. 18, 5321Ð5333.
[17] Burkhart, D., Sage, J., 2008. Cellular mechanisms of tumour suppression by the retinoblastoma gene. Nat Rev Cancer 8, 671–82.
[18] Chaves, M., Sontag, E., Albert, R., 2006. Methods of robustness analysis for Boolean models of gene control networks. IEEE Proceedings-
Systems Biology 153, 154–167. doi:10.1049/ip-syb.
[19] Chinnam, M., Goodrich, D., 2011. RB1, development, and cancer. Curr Top Dev Biol 94, 129–69.
[20] Conley, C., 1978. Isolated Invariant Sets and the Morse Index. Regional conference series in mathematics, R.I.:American Mathematical
Society.
[21] Crama, Y., , Hammer, P., 2011. Boolean Functions: Theory, Algorithms, and Applications. Cambridge University Press.
[22] Crawford-Kahrl, P., Cummins, B., Gedeon, T., 2019. Comparison of two combinatorial models of global network dynamics. SIAM J. Appl.
Dyn. Syst. 18, 418–457.
[23] Cross, F.R., Buchler, N.E., Skotheim, J.M., 2011. Evolution of networks and sequences in eukaryotic cell cycle control. Phil. Trans. R. Soc.
B 366, 3532–3544.
[24] Cummins, B., Gedeon, T., Harker, S., Mischaikow, K., 2017. Database of dynamic signatures generated by regulatory networks (ds-
grn), in: Koeppl, J.F.H. (Ed.), Computational Methods in Systems Biology - 2017. Springer. chapter 19, pp. 300–308. doi:10.1007/
978-3-319-67471-1.
[25] Cummins, B., Gedeon, T., Harker, S., Mischaikow, K., 2018. Model rejection and parameter reduction via time series. SIAM Journal on
Applied Dynamical Systems 17, 1589–1616.
[26] Cummins, B., Gedeon, T., Harker, S., Mischaikow, K., Mok, K., 2016. Combinatorial Representation of Parameter Space for Switching
Systems. SIAM Journal on Applied Dynamical Systems 15, 2176–2212.
[27] Cury, J., Monteiro, P., Chaouiya, C., 2019. Partial order on the set of boolean regulatory functions. preprint 2019 .
[28] Edwards, R., 2001. Chaos in neural and gene networks with hard switching. Diff. Eq. Dyn. Sys. , 187–220.
[29] Edwards, R., Glass, L., 2006. A calculus for relating the dynamics and structure of complex biological networks., in: Berry,R.S. and Jortner,J.
(eds) Adventures in Chemical Physics, Wiley, Hobeken, USA. pp. 151–178.
[30] Fages, F., Soliman, S., 2008. Formal cell biology in biocham, in: M. Bernardo, P. Degano, G. Zavattaro (Eds.), Formal Methods for Compu-
tational Systems Biology, Springer-Berlin-Heidelberg. pp. 54–80.
[31] Filippov, A., 1988. Differential equations with discontinuous right-hand sides. Kluwer Academic Publishers Group.
[32] Fisher, J., Henzinger, T., 2007. Executable cell biology. Nat. Biotechnol 25, 1239–1250.
[33] Gedeon, T., Cummins, B., Harker, S., Mischaikow, K., 2018. Identifying robust hysteresis in networks. PLoS Comput Bio 14, e1006121.
doi:https://doi.org/10.1371/journal.pcbi.1006121.
[34] Gedeon, T., Harker, S., Kokubu, H., Mischaikow, K., Oka, H., 2017. Global dynamics for steep sigmoidal nonlinearities in two dimensions.
Physica D 339, 18–38.
[35] Glass, L., Kauffman, S., 1972. Co-operative components, spatial localization and oscillatory cellular dynamics. Journal of Theoretical Biology
34, 219–37. URL: http://www.ncbi.nlm.nih.gov/pubmed/5015702.
[36] Glass, L., Kauffman, S., 1973. The logical analysis of continuous, non-linear biochemical control networks. Journal of Theoretical Biology
39, 103–29. URL: http://www.ncbi.nlm.nih.gov/pubmed/4741704.
[37] Glass, L., Pasternack, J.S., 1978. Stable oscillations in mathematical models of biological control systems. Journal of Mathematical Biology
T Gedeon: Preprint submitted to Elsevier Page 14 of 16
Dynamics of networks
6, 207–223. URL: http://link.springer.com/article/10.1007/BF02547797.
[38] Gonzalez, A., Naldi, A., Sanchez, L., Thieffry, D., , Chaouiya, C., 2006. Ginsim: a software suite for the qualitative modelling, simulation
and analysis of regulatory networks. BioSystems 84, 91–100.
[39] Guziolowski, C., Videla, S., Eduati, F., Thiele, S., Cokelaer, T., Siegel, A., JSaez-Rodriguez, J., 2013. Exhaustively characterizing feasible
logic models of a signaling network using answer set programming. Bioinformatics 29, 2320–2326.
[40] Harker, S., 2017. Dsgrn software. https://github.com/shaunharker/DSGRN. doi:10.5281/zenodo.805413.
[41] Hodgkin, A., Huxley, A., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerve. The
Journal of Physiology 117, 500–544.
[42] Huttinga, Z., Cummins, B., Gedeon, T., Mischaikow, K., 2018. Global dynamics for switching systems and their extensions by linear differ-
ential equations. Physica D: Nonlinear Phenomena 367, 19–37. doi:10.1016/j.physd.2017.11.003.
[43] Ironi, I., Panzeri, L., 2009. A computational framework for qualitative simulation of nonlinear dynamical models of gene-regulatory networks.
BMC Bioinformatics 10, S14.
[44] Ironi, L., Panzeri, L., Plahte, E., Simoncini, V., 2011. Dynamics of actively regulated gene networks. Physica D: Nonlinear Phenomena 240,
779–794. URL: http://linkinghub.elsevier.com/retrieve/pii/S016727891000360X, doi:10.1016/j.physd.2010.12.010.
[45] de Jong, H., 2002. Modeling and simulation of genetic regulatory systems: a literature review. J Comput Biol 9, 67–103. doi:10.1089/
10665270252833208.
[46] de Jong, H., Gouze, J., Hernandez, C., Page, M., Sari, T., Geiselmann, J., 2004. Qualitative simulation of genetic regulatory networks using
piecewise-linear models. Bull Math Biol 66, 301–40. doi:10.1016/j.bulm.2003.08.010.
[47] Lander, A., 2014. A calculus of purpose. PLoS Biology 2, 0712.
[48] Leone, G., DeGregori, J., Sears, R., Jakoi, L., Nevins, J.R., 1997. Myc and ras collaborate in inducing accumulation of active cyclin e/cdk2
and e2f. Nature 387, 422–6.
[49] Manning, A., Dyson, N., 2012. RB: mitotic implications of a tumour suppressor. Nat Rev Cancer 12, 220–6.
[50] Mestl, T., Plahte, E., Omholt, S., 1995. A mathematical framework for describing and analyzing gene regulatory networks. J. Theor. Biol.
176, 291–300.
[51] Murrugarra, D., Veliz-Cuba, A., Aguilar, B., Laubenbacher, R., 2016. Identification of control targets in boolean molecular network models
via computational algebra. BMC Systems Biology 10, 94.
[52] Pardee, A., 1974. A restriction point for control of normal animal cell proliferation. PNAS 71, 1286–1290.
[53] Pauleve, L., Richard, A., 2012. Static analysis of boolean networks based on interaction graphs: A survey. Electronic Notes in Theoretical
Computer Science , 93–104.
[54] Perez-Roger, I., Kim, S.H., Griffiths, B., Sewing, A., , Land, H., 1999. Cyclins d1 and d2 mediate myc-induced proliferation via sequestration
of p27(kip1) and p21(cip1). EMBO J 18, 5310Ð5320.
[55] Poret, A .and Guziolowski, C., 2018. Therapeutic target discovery using boolean network attractors: improvements of kali. R. Soc. open sci
5, 121852.
[56] Richard, A., Comet, J., Bernot, G., 2005. R. thomas’ modeling of biological regulatory networks: introduction of singular states in the
qualitative dynamics. Fundamenta Informacae XX, 1–18.
[57] Saadatpour, A., Reka, A., 2013. Boolean modeling of biological regulatory networks: A methodology tutorial. Methods 62, 3–12.
[58] Saari, D., 1977. A global existence theorem for the four-body problem of newtonian mechanics. J. Differential Equations 26, 80–111.
[59] Sears, R., Nevins, J., 2002. Signaling networks that link cell proliferation and cell fate. J Biol Chem 277, 11617–11620.
[60] Snoussi, H., 1989. Qualitative dynamics of piecewise-linear differential equations: a discrete mapping approach. Dynamics and Stability of
Systems 4, 565–583.
[61] Snoussi, H., Thomas, R., 1993. Qualitative dynamics of piecewise-linear differential equations: a discrete mapping approach. Bull. Math.
Biol 55, 973–991.
[62] Terfve, C.e.a., 2012. Cellnoptr: a flexible toolkit to train protein signaling networks to data using multiple logic formalisms. BMC Syst. Biol
6, 133. doi:https://doi.org/10.1186/1752-0509-6-133.
[63] Thieffry, D., Romero, D., 1999. The modularity of biological regulatory networks. BioiSystems 50, 49–59.
[64] Thieffry, D., Thomas, R., 1995. Dynamical behaviour of biological regulatory networks- ii. immunity control in bacteriophage lambda. Bull.
Math. Biol 57, 277–297.
[65] Thomas, R., 1973. Boolean formalization of genetic control circuits. Journal of Theoretical Biology 42, 563–585. doi:10.1016/
0022-5193(73)90247-6.
[66] Thomas, R., 1991. Regulatory networks seen as asynchronous automata: A logical description. Journal of Theoretical Biology 153, 1–23.
doi:10.1016/S0022-5193(05)80350-9.
[67] Thomas, R., D’Ari, R., 1990. Biological Feedback. CRC-Press.
[68] Thomas, R., Thieffry, D., Kaufman, M., 1995. Dynamical behaviour of biological regulatory networks-i. biological role of feedback loops
and practical use of the concept of the loop-characteristic state. Bull Math Biol 57, 247–76.
[69] Tournier, L., Chaves, M., 2009. Uncovering operational interactions in genetic networks using asynchronous boolean dynamics. J Theor Biol
260, 196–209. doi:10.1016/j.jtbi.2009.06.006.
[70] Tournier, L., Goelzer, A., Fromion, V., 2017. Optimal resource allocation enables mathematical exploration of microbial metabolic configu-
rations. J. Math. Biol. 75, 1349–1380.
[71] Veflingstad, S.R., Plahte, E., 2007. Analysis of gene regulatory network models with graded and binary transcriptional responses. Biosystems
90, 323–339.
[72] Veliz-Cuba, A., 2012. An algebraic approach to reverse engineering finite dynamical systems arising from biology. SIAM J. Dynam. Sys. 11,
31–48.
[73] Videla, S., Konokotina, I., Alexopoulos, L., Saez-Rodriguez, J., Schaub, T., Siegel, A., Guziolowski, C., 2015. Designing experiments to
T Gedeon: Preprint submitted to Elsevier Page 15 of 16
Dynamics of networks
discriminate families of logic models. Frontiers in Bioengineering and Biotechnology 3, 131.
[74] Xia, Z., 1992. The existence of noncollision singularities in newtonian systems. Annals of Mathematics 135, 411–468.
[75] Xin, Y., Cummins, B., Gedeon, T., 2019. Multi-stability in epithelial-mesenchymal transition network. preprint .
[76] Yao, G., Lee, T., Mori, S., Nevins, J., You, L., 2008. A bistable Rb-E2F switch underlies the restriction point. Nat. Cell Biol 10, 476–482.
[77] Yao, G., Tan, C., West, M., Nevins, J., You, L., 2011. Origin of bistability underlying mammalian cell cycle entry. Molecular Systems Biology
485.
[78] Yousefi, M., Datta, A., Dougherty, E., 2012. Optimal intervention strategies for therapeutic methods with fixed-length duration of drug
effectiveness. Signal Process IEEE Trans 60, 4930–44.
[79] Yousefi, M., Datta, A., Dougherty, E., 2013. Optimal intervention in Markovian gene regulatory networks with random-length therapeutic
response to antitumor drug. Biomed Eng IEEE Transac 60, 3542–52.
[80] Yousefi, M., Dougherty, E., 2013. Intervention in gene regulatory networks with maximal phenotype alteration. Bioinformatics 29, 1758–67.
[81] Yousefi, M., Dougherty, E., 2014. A comparison study of optimal and suboptimal intervention policies for gene regulatory networks in the
presence of uncertainty. EURASIP J Bioinforma Syst Biol 2014, 6–6.
T Gedeon: Preprint submitted to Elsevier Page 16 of 16