The goal of this paper is to extend an existing dynamical systems method for modeling transcriptional regulation in gene networks to include protein transport regulation. The framework under study is a combined discrete and continuous approach implemented as a numerical method called DSGRN (Dynamic Signatures Generated by Regulatory Networks) [9] and we will refer to both the theoretical framework and the software as DSGRN. DSGRN is a modeling technique that uses combinatorial dynamics to fully describe the range of long-term dynamics that is compatible with a given network structure across all parameters [9,8,23]. This approach has been used successfully to examine the robustness of different phenotypes across parameter space in the cell cycle [24], the epithelial-mesenchymal transition in cancer [41], and synthetic biology design [20]. Motivated by switching system network models [38,39,12,16,30,31], which use piecewise constant nonlinearities and are themselves motivated by Boolean models [3], DSGRN was originally designed to model transcriptional interactions in biological networks. In these interactions, the rate of production of one quantity (y) is regulated by the presence of another quantity (x). However, there are other types of interactions that are modeled by regulatory networks, in particular, when (1) the activity rather than the quantity of y is regulated by x, (2) the decay rate of y is regulated by x, and (3) the transport of y between compartments is regulated by x. Regulation of activity and decay rate has been recently incorporated into DSGRN [11]; here we seek to further extend the methodology to transport regulation.
A fundamental assumption of DSGRN is that an edge from node x to node y in a regulatory network represents two levels of influence, L y,x < U y,x . When the concentration of x is higher/lower than a threshold ζ y,x , one of these two levels of influence affects y. The position of x with respect to ζ y,x is a combination of levels of its own inputs, where the amount of x increases if the input combination exceeds the decay rate γ x of x, and decreases otherwise. Therefore an arbitrary edge x to y has three associated parameters ζ y,x , L y,x , U y,x . In addition, each node x has an associated decay parameter γ x . The union of these non-negative parameters forms a parameter space.
Evolution of concentrations can be coarsely represented by transitions in a state transition graph (STG). The thresholds ζ y,x naturally decompose phase space of concentrations into rectangular domains. Each domain is represented by a node in the STG, i.e. a state. For each parameter set, transitions in the STG match increasing/decreasing behavior for each x in the particular state. An STG can be thought of as a multi-valued map on the set of states, or, equivalently, as a graph with nodes corresponding to states. The long-term behavior of the dynamics of an STG is captured by computing its graph of strongly connected path components, ordered by reachability, called the Morse graph. The key observation is that STGs, and hence Morse graphs, only depend on a finite collection of potential inequalities between the input levels to a node x and the collection of its output node thresholds. Since all such collections of inequalities can be enumerated, one can decompose parameter space into a finite number of regions called DSGRN parameter nodes and capture the entire collection of possible Morse graphs for a given network. This leads to a comprehensive (though coarse) description of all dynamical behavior that a network structure can support through different parameterizations.
The DSGRN approach is closely related to the idea of logical bifurcation analysis introduced in [1]. In particular, the parameter graph captures proximity relationships of domains in parameter space and each node of the parameter graph can be viewed as a collection of monotone Boolean maps [7]. Therefore any path in the parameter graph is a one-dimensional logical bifurcation diagram of [1].
There are three challenges associated to the extension of DSGRN to include transport terms. The first is the loss of a product structure of the parameter graph. In particular, the DSGRN parameter graph P G = Π n i=1 P G(x i ) is a product of factor graphs P G(x i ), each of which captures the relationship between input levels of influence and output thresholds at a single node x i . Transport regulation always affects two compartments/nodes, as a decrease of the transported term in one compartment is balanced by an increase in the target compartment. This leads to loss of product structure in the parameter graph between the factor graphs corresponding to these two nodes through the introduction of dependencies. This requires a fundamental modification to the idea of the DSGRN parameter graph. The second challenge is that DSGRN is well-defined only if all thresholds are distinct. This means that if a variable corresponding to node x crosses a threshold, this may affect the dynamics of at most one other node y. However, when a transport regulation variable x crosses a threshold, it affects two variables simultaneously. These variables represent the concentration of the transported molecule in the compartment from which the transport originates and the concentration in the compartment where the transport terminates. These first two challenges also arise in more general biochemical networks where enzymatic reactions convert a substrate to a product. Then the decrease of substrate concentration is accompanied by the proportional increase in the product concentration and this process can be thought of as a generalized transport between these two compartments. Lastly, DSGRN was originally formulated for graphs without multi-edges, so that only one edge can exist from a source to a target. While not strictly a necessary consequence of the presence of a transport term, a protein x can act as a transcriptional regulator of y and at the same time, control transport of y between compartments. We will see this in our example system, which leads us to consider the possibility of two control edges between x and y. In this paper we will show how to address these three issues and extend DSGRN to transport regulation interactions.
We illustrate our approach by analyzing a p53-Mdm2 three node network model, which exhibits a mixture of transcriptional regulation, activity regulation and a prominent regulated transport term. Protein p53 is a transcription factor that responds to cell damage by regulating genes involved in growth arrest, DNA repair, apoptosis and cellular senescence [21,35,40]. Mdm2 is a ubiquitin ligase that increases the degradation rate of p53 when in the nucleus. There are two influences of p53 on Mdm2, a transcriptional up-regulation and transport down-regulation from the cytoplasm to the nucleus. Nuclear and cytoplasmic Mdm2 together with p53 then form a simplified threenode negative feedback network (see Figure 3b and Figure 3c) that models the control of the cell's response to damage. This network was analyzed by a switching system model in [2], which followed an earlier work using traditional ODE models [1,36]. The goal of these papers was to reproduce two types of oscillations observed experimentally in [25] under two different conditions: low and high irradiation of the cell, intended to produce lower and higher levels of cell damage. It was found that the two radiation levels were associated to two different frequencies of response, lower (higher) frequency for lower (higher) cell damage. In [2], a combination of switching system modeling and extensive numerical simulations demonstrated that these oscillations can coexist as stable periodic orbits in a particular parameter regime.
We validate our extension of DSGRN to transport regulation by replicating these results. We show the coexistence of higher and lower frequency oscillations within the state transition graph, a structure that we call A 8 because it resembles a figure 8, for some DSGRN parameter nodes. We then sample these nodes, simulate ODE models with Hill function nonlinearities with sampled parameters, and find both types of stable periodic orbits at nearby, although not identical, parameter sets. While we do not find real-valued parameter samples where both periodic orbits coexist, this is not surprising given that the parameter region of coexistence found in [2] is very narrow. Moreover, the coexistence of two orbits implies that a change in initial conditions, i.e. protein concentrations, is enough to trigger a shift between orbits. However, we note that the distinct periodic cycles seen in the data induced by low and high radiation levels may be realistically modeled as a parameter change rather than as a change in initial conditions. For these reasons, we view the existence of both periodic orbits for parameter values that are near each other as a successful reproduction of the experimental data.
The theoretical work presented here is intended to be incorporated into the DSGRN software in the future. This will automate the manual work in Section 4 and make this type of transport modeling accessible to other researchers. Since DSGRN calculations of the state transition graphs for the entire parameter space are almost instantaneous for small network models, DSGRN is a useful tool for quickly scanning the parameter space of a network model to locate complex dynamic signatures such as two stable periodic orbits with particular shapes. This allows the comprehensive assessment of the validity of hypothesized network models after experimentation, as well as prediction of expected behavior before experimentation.
DSGRN is built on the ideas of Boolean modeling of gene regulatory networks that originated with Thomas, Kauffman and Glass [26,27,38,38]. To make the discussion in following sections more accessible, we illustrate the development of these ideas on a simple two node network example in Figure 1a. To such a network, following [26,38,12,13,16,30] one associates a system of ODEs
These nonlinearities are piecewise constant. The values of L < U and ζ are parameters that can be easily translated into parameters of a family of Hill function models via
parameterized by the Hill coefficient n. In what follows, in order to distinguish functions assigned to different edges, we will use the subscript (2, 1) to denote the function σ + 2,1 associated with edge 1 → 2, as well as all corresponding parameters U 2,1 , L 2,1 , ζ 2,1 .
The dynamics of the ODE model of the network (1) are relatively easy to analyze. The thresholds ζ 1,1 , ζ 2,1 separate the positive y 1 axis into three intervals and ζ 1,2 separates the positive y 2 axis into two intervals. Therefore the positive orthant R 2+ is divided into 6 domains (Figure 1b). In each of them the dynamics are simple; the flow lines are linear and converge toward a target point, the values of which are shown in each domain in Figure 1b. We want to focus on two issues: capturing the dynamics in the form of a state transition graph (STG) and the dependence of such dynamics on parameters of the system. The normal components of the vector field across domain boundaries can be analytically calculated from the structure of (1), which captures the flow from domain to domain, see Figure 1c. The vector field will flow from one domain to a neighboring domain if the domain and its target point differ in that coordinate. This represents asynchronous update dynamics. To define the STG we represent each of the 6 domains in R 2+ by a state in the state space X = {0, 1, 2} × {0, 1}; the integer in each component denotes the interval along the corresponding axis. The edges of the STG represent all possible transitions induced by the dynamics to neighboring domains, see Figure 1d. The parameters of system (1) are two decay parameters that can be associated with the nodes of the network, γ 1 , γ 2 , and three parameters L, U, ζ for each of the three edges. Therefore the parameter space for system (1) is the positive orthant R 11+ . Since the thresholds Th(1) := {ζ 1,1 , ζ 2,1 } are both thresholds for variable y 1 , there are two possible orders,
} be the set of edges that start at node 1, where e i,j is an edge (j → i) from node j to node i. Note that since each threshold is associated with exactly one edge, there is a bijection between T(1) and Th (1). With an eye toward notation used in the extension of DSGRN in Section 3, we describe each of the two orders of the set Th(1) by maps θ 1 : {1, 2} → T(1) and φ 1 : {1, 2} → T(1), corresponding to the orders ζ 1,1 < ζ 2,1 and ζ 2,1 < ζ 1,1 respectively. The maps are given by ( 2)
The order θ 2 : {1} → T(2) trivially consists of one element.
Letting φ 2 := θ 2 , we call the collections of maps θ := (θ 1 , θ 2 ) and φ := (φ 1 , φ 2 ) order parameters.
Note that the STG dynamics are directly derived from knowledge of the target point of each domain or, more precisely, knowledge of the domain in which the target point lies. The target point of the solution to (1) in each particular domain is
evaluated at any point (y 1 , y 2 ) in the domain. Selecting for the sake of an example the first order θ, Figure 1b depicts the target point values 1 γ 1 σ - 1,2 (y 2 )σ + 1,1 (y 1 ) and 1 γ 2 σ + 2,1 (y 1 ), where the subscripts on L, U have been suppressed for clarity. Observe that the identity of the domain of the target 1 2 (0, 0)
Figure 1 point depends on two sets of inequalities: in the second component either
Since there are only two possible values L 2,1 < U 2,1 of σ + 2,1 (y 1 ), there are three regions in the parameter space
Notice the swap in subscript order between ζ 1,2 and L 2,1 , U 2,1 . This is not an error; the target point constraints depend on the ordering of the input values of y 2 (that is L 2,1 , U 2,1 ) to the output threshold ζ 1,2 of y 2 .
We encode the inequalities via collections of monotone Boolean functions, one for each out-edge of a node. In the case of y 2 , there is a single out-edge 2 1 and we define the associated monotone Boolean function f 1,2 : I 2 → B, where I 2 is a Boolean representation of the set of admissible input words {L 2,1 , U 2,1 } via bijection S(L 2,1 ) = 0 and S(U 2,1 ) = 1. The Boolean value of f 1,2 indicates whether that word lies below or above the threshold. In the range off 1,2 we set f 1,2 (S(q)) = 0 when q < γ 2 ζ 1,2 and f 1,2 (S(q)) = 1 when q > γ 2 ζ 1,2 , for q ∈ {L 2,1 , U 2,1 }. Therefore the three choices above correspond to three Boolean functions
Each of f 2 := (f 1,2 ), g 2 := (g 1,2 ), and h 2 := (h 1,2 ) is a logical parameter at node 2.
Similar considerations about the position of the first coordinate of the target point leads to a set of logical parameters that describe all possible positions of partially ordered values
within three intervals (0, 1b, the flow direction in Figure 1c, and the STG in Figure 1d correspond to a choice of parameters for the differential equation (1) which satisfies
The Boolean function representation for node 1 for this choice is given by
Each of 20 distinct possible choices for the function
) is a logic parameter f 1 (•; θ) which depends on the order of thresholds and thus the order parameter θ. The reduction of 64 potential pairs of Boolean functions consistent with positive regulation down to 20 happens due to a requirement of monotonicity and a deeper idea of realizability discussed in [9].
The collection P = {(θ, (f 1 , f 2 ))} of paired order and logic parameters describes 120 = 2 * 20 * 3 disjoint regions of decomposition of the parameter space R 11+ , each of which determines an STG. We organize this information in the form of a parameter graph (P G) whose nodes are P and edges represent two regions that differ by a single inequality. We provide a more precise definition in the next section. The parameter graph P G indexes information about the entire repertoire of dynamics of the network as described by STGs. In particular, this structure can be searched for prevalence of certain dynamic behaviors (i.e. types of attractors) or changes in dynamics as a particular parameter (or a set of parameters) changes. The work [1] used the term logical bifurcation diagram to describe changes in the STG as a function of changing parameters. Any logical bifurcation diagram can be constructed from the parameter graph by a projection.
We close this section with a few remarks. The size of the parameter graph grows very quickly with the size of the network, in particular, the number of edges. However, note that logic parameters f 2 and h 2 shown in (3) impart trivial functionality to edges in the network. Regardless of the input from y 1 , under f 2 (h 2 ) the variable y 2 can only decrease (increase) through the threshold ζ 1,2 ; i.e. node 2 is non-responsive to node 1 at logic parameters f 2 , h 2 . Furthermore, y 2 can never nontrivially impact the dynamics of y 1 , also because it only has unidirectional motion across the threshold ζ 1,2 ; i.e. 2 is an ineffective regulator of 1 at f 2 , h 2 . On the other hand, the choice of g 2 in our example represents a situation where 2 is responsive to 1 and is also an effective regulator of 1. When node i is responsive to all of its inputs and is an effective regulator of all of its outputs at some logic parameter f i , we say that f i is an essential parameter for i or simply essential. Given a parameter p = (θ, f ) ∈ P, if every component f i of the logic parameter f is essential, then we say that p is an essential parameter for the network or simply essential parameter. This is explained in more detail in Section 3.5.
The second comment is on whether STG is always well defined for every p ∈ P. The only problem that may occur is when the network has a negative self-edge on some node i [9]. In this case, in two domains bordering ζ i,i the vector field can point toward ζ i,i , which signals a potential problem with forward continuation of (1) or at least a presence of flow along the boundary instead of transverse flow across the boundary [13]. In this situation the border along ζ i,i between two domains is called a black wall [37,16]. An important observation is that under a generic assumption that all thresholds ζ j,i are distinct for all target nodes j a of a node i, the black walls are consistent: all other components besides y i do not change on ζ i,i . Therefore one can define a consistent vector field along a consistent black wall.
Unfortunately, the presence of the transport regulation induces an identical term in two equations that describe donor and recipient compartments of the transport. This means that any threshold in the transport regulation term will affect both equations. Therefore two components of the vector field may change on the same threshold, which may lead to inconsistent black walls. We describe a method in Section 3.7 that extracts a collection of valid STGs in this situation.
In this section we describe a construction of the DSGRN parameter graph and dynamics for transport networks that incorporates both regulatory effects between nodes and transport between compartments. Unlike traditional DSGRN regulatory networks, transport networks allow for (1) multiple edges between nodes and (2) pairs of edges associated with transport or regulation of transport that have the same threshold. Transport networks are formally described in Section 3.1.
The order parameter, logic parameter, parameter graph, and essential parameter graph for a transport network are described in Sections 3.2, 3.3, 3.4, and 3.5, respectively. These are defined differently in the transport extension of DSGRN, due to the fact that two out-edges of a single node i have identical thresholds, ζ j,i = ζ i,i , one of which is always associated to a black wall (Section 3.6). This gives rise to the aforementioned inconsistent black walls. We resolve the inconsistency by perturbing the two thresholds so that they are not equal: either ζ i,i < ζ j,i or ζ j,i < ζ i,i . This results in two unfolded parameters for each unperturbed parameter, one for each ordering (Section 3.7).
The most important notation for this section is established in Table 1. Note that for a network edge e given by i → j, the subscript convention (j, i) will sometimes be replaced by e, especially when the nodes in the edge are not explicitly identified. As an example, instead of writing σ + j,i , we may instead write σ + e . This convention is used in Table 1. 3.1. The transport network. We consider here actively regulated, unidirectional transport of a substrate y from a compartment c i to a compartment c j . We denote the concentration of y in compartment c i as y i and similarly y j in c j . We model this situation by a pair of differential equations of the form
where Λ i and Γ ji are combinations of the discontinuous switching functions σ + and σ -introduced in Section 2. Note that for any i the term Λ i or Γ ji may be zero, but if Γ ji is not zero, it must appear in both i-th and j-th equation. The parameter A > 0 reflects the difference in the volumes of compartment c i and c j . The multiplier ± σ + (y i ) in the transport termsσ + (y i )Γ ji and + Aσ + (y i )Γ ji indicate that as the concentration y i increases, the transport of y across the
set of |S| dimensional Boolean vectors indexed by the elements of the set S TN = (V, E, E tr ) regulatory network with nodes V , edges E, and transport edge pairs E tr i, j, u, v ∈ V nodes in the network with integer values 1, . . . , N y i , y j ∈ R + variables representing concentrations in the switching system for nodes i, j T(v) set of target edges from node j T r (v) the subset of T(v) that are not transport edges T tr (v) the set of transport edge pairs involving two edges in
the subset of S(v) that are not transport edges S tr (v) the set of transport edge pairs such that one of the edges is in
the lower and upper real-valued response parameters associated to an edge e ∈ E 0 < ζ e a real-valued threshold associated to edge e ∈ E θ v a map that orders the thresholds ζ e for e ∈ T(v) θ an order parameter, (θ
the set of admissible words, v) , and the domain of
DSGRN parameter graph with nodes P X integer state space of the DSGRN switching system x ∈ X a possible integer state of the nodes in TN, x = (x 1 , . . . , x N ) F(x; p) a multi-valued map on X that describes the dynamics of TN at a parameter p ∈ P Table 1. Selected Notation membrane also increases. Note that the sets of nodes {i 1 , . . . , i t } that regulate abundance and activity of y i , the set of nodes {j 1 , . . . , j v } that regulate abundance and activity of y j , and the set of nodes {u 1 , . . . , u s } that regulate the transport may have non-empty intersections.
We define a transport network TN = (V, E, E tr ) as consisting of nodes V , signed directed edges E ⊂ V × V , and unordered transport edge pairs E tr ⊂ E × E. With regard to (5), the TN transport nodes i, j ∈ V are a transport pair, the nodes {u 1 , . . . , u s } are transport regulators, and the nodes {i 1 , . . . , i t } and {j 1 , . . . , j v } are transcriptional regulators, where a node k may be both a transcriptional and a transport regulator. We denote the sign of an edge e ∈ E by sgn(e) ∈ {-1, 1}. Each edge e ∈ E corresponds to an instance of a switching function σ ± e in an ODE system representing TN that contains at least one pair of equations taking the form of ( 5). The superscripts reflect the sign of the edge e = (u, v), with the exception of the expression -σ + e , where sgn(e) = -1, and -σ - e , where sgn(e) = 1. The edge sign therefore captures whether increasing the concentration of the source node u tends to increase (sgn(e) = 1) or decrease (sgn(e) = -1) the concentration of the target node v. Since the transcriptional and transport regulators may have non-empty intersection, E may include multi-edges. The purpose of the transport edge pairs E tr is to indicate when the same switching function σ ± (y u k ) or σ ± (y i ) appears in both equations for ẏi and ẏj . That is, the pairs of transport edges E tr correspond to the transport terms in (5). The transport edges satisfy the following properties (1) Every pair consists of unique edges: if {e 1 , e 2 } ∈ E tr then {e 1 , e} / ∈ E tr for all e ∈ E \ {e 2 } and {e 2 , e} / ∈ E tr for all e ∈ E \ {e 1 }.
(2) If {e 1 , e 2 } ∈ E tr then e 1 = (u, v 1 ) and e 2 = (u, v 2 ) have the same source node u.
(3) If {e 1 , e 2 } ∈ E tr , then e 1 and e 2 have opposite sign: sgn(e 1 ) = -sgn(e 2 ). This reflects the fact that -σ + (y i )Γ appears in the equation for ẏi but +σ + (y i )Γ j,i appears in the equation for ẏj in (5).
Properties ( 1) and ( 2) will be essential to describe DSGRN parameters and dynamics. Property (3) is a property of the model, but is not a necessary assumption to produce dynamics.
Notice that this formulation permits the existence of multi-edges. For example, in (5) there may be two edges from node u k to i, one of which regulates transcription while the other regulates transport.
3.2. The Order parameter. Let T(u) denote the set of target edges that begin at node u.
In traditional DSGRN (networks without transport pairs where E tr = ∅), the order parameter θ = (θ 1 , . . . , θ N ) is a collection of invertible maps where the map θ v for node v describes an ordering of T(v), see ( 2) given in the introductory example. This ordering reflects the ordering of the realvalued thresholds Th(v) := {ζ e } e∈T(v) which appear in the switching functions σ ± e corresponding to each edge. For the transport network TN, the number of unique values in Th(v) may not be the same as the number of edges in T(v) because transport edge pairs have the same threshold by definition. An order parameter for a transport node is defined to reflect this. Let T r (v) ⊂ T(v) denote the set of target edges of v which do not appear as part of a transport edge pair in E tr . Let T tr (v) denote the transport edge pairs {e 1 , e 2 } with e 1 ∈ T(v). Note that the second property of transport edges states that e 1 ∈ T(v) implies e 2 ∈ T(v). The order parameter for the transport network is a collection θ = (θ 1 , . . . , θ N ) where for each node v, θ v is a choice of bijection
. We let Θ denote the set of all order parameters for TN.
| in the presence of transport regulation by property (1) of transport networks. In particular, each edge in T(v) \ T r (v) occurs in exactly one transport edge pair, so two edges are counted as one.
In the original DSGRN formulation the order parameter for node u has a domain of {1, . . . , |T(v)|}. Due to the change in the domain size of the transport order parameter, we also redefine the state space of the DSGRN switching system, previously introduced via example in Section 2. The original DSGRN definition of the state space is X = v∈V {0, 1, . . . , |T(v)|}, and the new definition is
The state space both labels the nodes of the STG and provides a signature of the domain in phase space where a point in R N + is located with respect to the thresholds {ζ e } e∈E . We call an element x ∈ X a state.
The logic parameter. Recall from Section 2 that a logic parameter at a network node v is a collection of monotone Boolean functions (MBFs), one for each target edge of v. In the original DSGRN formulation, the domain of the MBFs is always B S(v) , the set of |S(v)| dimensional Boolean vectors with components indexed by the edges in S(v). This section will carefully build a restricted domain I v (θ) ⊂ B S(v) for the MBFs of the transport logic parameter that is dependent on the choice of order parameter, and then formally introduce the resulting MBFs.
Let S(v) denote the set of source edges that end at node u, S r (v) ⊂ S(v) denote the set of sources edges of v which do not appear as part of a transport edge pair in E tr , and S tr (v) denote the transport edge pairs {e 1 , e 2 } such that at least one edge is a source edge of j, i.e. either e 1 ∈ S(v) or e 2 ∈ S(v). The effect of the source edges S(v) on v is described by an input word that is a Boolean vector b ∈ B S(v) . These input words are arguments for the logic parameter f v associated to node v, which was informally introduced in Section 2 and will be formally defined below. The logic parameters encode which inputs tend to move y v above or below the thresholds associated to output edges of node v in solutions to the ODE system containing (5).
To describe the input word associated with a state x ∈ X, fix an order parameter θ ∈ Θ and consider an edge e = (u, v) ∈ S r (v). The input from node u to node v via e is defined by b e (x; θ) := 0, x u < θ -1 u (e) and sgn(e) = 1 or θ -1 u (e) ≥ x u and sgn(e) = -1 1, θ -1 u (e) ≥ x u and sgn(e) = 1 or x u < θ -1 u (e) and sgn(e) = -1
.
In other words, we compare the state x u to the value of θ -1 u (e), which represents the ordered position of the threshold ζ e along the y u axis in phase space. Depending on the sign of the edge, the binary output of b e (x; θ) indicates whether the input to node v along the edge e is low (b e = 0) or high (b e = 1) at each state.
For an edge e = (u, v) ∈ S(v) \ S r (v) which is part of a transport edge pair, let {e, e } ∈ S tr (v) be the pair. The input from node u to node v via this transport edge e is defined in the same way as above. The only difference is that the edge e is not in the range of θ u , but only the pair {e, e } is. The input is given by b e (x; θ) := 0, x u < θ -1 u ({e, e }) and sgn(e) = 1 or θ -1 u ({e, e }) ≥ x u and sgn(e) = -1 1, θ -1 u ({e, e }) ≥ x u and sgn(e) = 1 or x u < θ -1 u ({e, e }) and sgn(e) = -1
.
The input due to all source edges to node v at x ∈ X is then collected to form the input word which is given by the Boolean vector
The set of admissible words for node v and order parameter θ are
Remark 3.1. Which words belong to the admissible set I v depends on the constraint arising from the existence of multiple edges between node u and node v. For illustration, assume that there are two edges {e 1 , e 2 } = S(v) both of which initiate at the same node u and terminate at v. Assume that the corresponding thresholds
). Then there is no y ∈ R + that simultaneously satisfies
This makes one of the four words {00, 01, 10, 11} = B S(v) not admissible; for example, if both edges are positive then the word b = b e 1 b e 2 = 01 is not admissible. When all edges e ∈ S(j) have distinct initial nodes, there are no constraints and
as in traditional DSGRN and, in fact, I v (θ) = I v does not depend on the order parameter θ.
Let {e 1 , e 2 } ∈ E tr be a transport edge pair. By property (3) of transport edges, b e 1 (x; θ) = ¬b e 2 (x; θ) for all x since e 1 and e 2 have opposite signs. Taking into account this constraint, the set of admissible sentences I(θ) is the collection of admissible words
e 2 for all (e 1 = (u, v 1 ), e 2 = (u, v 2 )) ∈ E tr . This completes the construction of the domain I(θ) of the monotone Boolean functions that comprise a logic parameter. The construction of the order-parameter dependent MBFs themselves now follows.
We will first define an MBF f q for each edge or transport edge pair q ∈ T r (j) ∪ T tr (j) for a node j. An MBF for the node v, f v , is composed of these f q , and finally the collection of f v over all nodes v ∈ V forms the logic parameter f . Definition 3.2. Let ≺ denote the partial order of the Boolean lattice B n and let I ⊂ B n be a sublattice. A function g :
In particular, note that I v (θ) ⊂ B S(v) is a sublattice with partial order ≺ inherited from B S(j) .
For each edge or transport edge pair, q ∈ T(v) r ∪ T tr (v), originating in node v, we let f q (•; θ) : I v (θ) → B be a positive monotone Boolean function such that
for all b ∈ I v (θ). Define the collection of positive MBFs corresponding to edges whose source node is v:
Finally, define the aggregation of all MBFs over the nodes of the network:
For clarity and brevity we will use the notation f (θ) to denote this collection. We call f (θ) a logic parameter for order parameter θ and let L(θ) be the set of all such logic parameters.
3.4. The parameter graph. The set of parameter nodes is given by the set of pairs of order parameters and logic parameters:
Two parameter nodes (θ, f (θ)), (φ, g(φ)) ∈ P are adjacent when exactly one of the following conditions is satisfied.
• Order adjacency: The order parameters θ and φ differ only in one pair of neighboring entries. Explicitly, there is a unique node v, and an integer k ∈ {1, . . .
In addition, the maps f (θ) and g(φ) are the same in the following sense.
(1) When the edges or transport edge pairs θ v (k) and θ v (k + 1) have different target nodes for every edge, the admissible sentences I(θ) and I(φ) are identical by Remark 3.1. Therefore f (θ) and g(φ) have the same domain. In this case we require f (•; θ) = g(•; φ).
(2) If θ v (k) and θ v (k + 1) have a common target w, then I w (θ) = I w (φ). In this case we require f u (•; θ) = g u (•; φ) for u = w and f w (b w (x; θ); θ) = g w (b w (x; φ); φ) for all x ∈ X.
(3) Finally, we require that the logic functions for the edges or transport edge pairs θ v (k) and θ v (k + 1) agree:
which induces g φv(k) = g φ j (k+1) .
• Logical adjacency: θ = φ and the logic parameters f and g differ at a single input. Explicitly, there is a unique node v, unique q ∈ T r (v)∪T tr (v), and unique b ∈ I v (θ) = I v (φ) such that f q (b) = g q (b). For an edge or transport edge pair q = q or input b = b, we require
The parameter graph is the undirected graph P G = (P, E) whose nodes are parameter nodes and edges are given by adjacency.
We remark that the set of parameter nodes P for a transport network as defined here avoids the notion of realizability as introduced in [9]. A parameter node p ∈ P is realizable if, given the algebraic structure of the nonlinearities Λ i and Γ ji in the ODE system underlying TN, the parameters L, U , ζ, and γ can be chosen to be consistent with the order and logic parameters of p. This is a deep problem in algebraic geometry that so far has only been computationally solved for limited cases in DSGRN [32]. In our formulation, P may contain unrealizable parameter nodes.
A parameter graph contains nodes p = (θ, f (θ)) where f q (•; θ) : I v (θ) → B is constant for some node v and some edge or transport edge pair q. In this situation, the edge(s) q is (are) either always "active" when f q (•; θ) ≡ 1 or always "inactive" when f q (•; θ) ≡ 0. This means that the dynamics of the network, as we will define in Section 3.6, are not dependent on this edge. Therefore, P G represents not only the dynamics of a network, but also the dynamics of all proper subnetworks resulting from removing edges which do not affect their target nodes.
The subgraph of P G that represents the dynamics of the full network is called an essential parameter graph and its nodes essential parameter nodes. The parameter node p = (θ, f (θ)) is output essential if for each edge or transport edge pair q ∈ T r (v) ∪ T tr (v) emanating from node v, the function f q (θ) : I v → B is non-constant. The parameter node p = (θ, f (θ)) is input essential if for each node v and edge e ∈ S(v), there is a q ∈ T r (v) ∪ T tr (v) and words a, b ∈ I v (θ) satisfying a e = b e but a = b for all = e such that f q (a) = f q (b).
The parameter node p = (θ, f (θ)) is essential if it is both input essential and output essential.
3.6. DSGRN dynamics. We now describe a family of multi-valued maps F : X × P ⇒ X where recall that the state space X is given by
These parameter-dependent multi-valued maps capture the discrete dynamics of the DSGRN switching system for a transport network TN = (V, E, E tr ). They are a generalization of the state transition graphs introduced in Section 2.
Let (θ, f (θ)) ∈ P be given and consider an admissible sentence (b 1 , . . . , b N ) ∈ I(θ), where each b v (x; θ) is an admissible input word to the node v realized at state x ∈ X.
Definition 3.3. Let (θ, f (θ)) ∈ P.
(1) The synchronous Boolean dynamics F 0 : X × P → X is defined by
The state F 0 (x; (θ, f (θ)) is the target for state x and is analogous to the target point for the ODE formulation described in Section 2.
(2) The asynchronous Boolean dynamics F : X × P ⇒ X is defined by
• For any v and
DSGRN is based on asynchronous Boolean dynamics. The multivalued map F : X × {p} ⇒ X is often represented as a state transition graph (STG) with nodes (states) X and edges E X , where an oriented edge x → x ∈ E X exists if and only if x ∈ F(x; p). These edges only exist between states that differ by 1 in a single direction. Each state x corresponds to a domain in the phase space R n bounded on all sides by thresholds, zero, or extending to infinity. The edges in an STG represent the direction of transition between two domains consistent with local solutions of a switching ODE [10,27,16], as described in (1).
Consider two neighboring states x, x that differ in u-th direction with x u = x u + 1. When there is a unique threshold ζ e , associated to an edge e = (u, v) ∈ E in the transport network, between domains corresponding to x and x , only the v component F v (x) may differ from F v (x ). In particular, F u (x) = F u (x ) =: k. If this common value k > x u then STG contains x → x ∈ E X and if k < x u , then the STG contains x → x ∈ E X . It follows that for u = v, at most one of the edges x → x and x → x is present in the STG.
However, if the edge e = (u, u) ∈ E is a negative self-loop, both edges x → x and x → x may be present, see Figure 2. We call this situation a black wall, which refers to an image of a hyperplane corresponding to y u = ζ u,u in the phase space which is approached from both sides by solutions of switching ODEs. Black walls present technical challenges for the existence and uniqueness of solutions of the switching ODE [18,13], in some cases this inconsistency can be resolved. There is considerable literature on regularization both in more general systems [?] and in the switching systems considered here. For the switching systems multiple papers use limits of trajectories of the smooth perturbed systems to define solutions of the switching system [30, 31, ?, ?]. In our approach, we distinguish two types of black walls. Suppose θ u (x u ) = e ∈ T r (u) and there is no transport edge pair of the form {(u, u), e } ∈ T tr (u), so that the domains corresponding to x and x are separated by a single threshold. This case is called a consistent black wall (Figure 2a) and in this case an STG can be defined as follows. When x u = k, we introduce a new state k * into the u-th component of X to form
It can be shown [20] that since only the u-th component of F is changing between x and x the STG can be extended to a new extended space X. On the other hand, the definition of an STG in case of inconsistent black walls (Figure 2b), which occur when θ u (x u ) = {e, e } ∈ T tr (u) is a transport edge pair, is more complicated.
Inconsistent black walls and parameter unfolding. The presence of the transport term always produces DSGRN parameters with an inconsistent black wall since transport from compartment c i to c j corresponds to transport edge pairs {e 1 , e 2 } = {(i, i), (i, j)} where the self-edge is negative: sgn(e 1 ) = -1. This can be seen in ( 5), by examining the terms ±σ + (y i ) in the equations for ẏj and ẏi , respectively. This produces adjacent states x, x ∈ X, which differ in only the ith component, such that x ∈ F(x ), x ∈ F(x), but F j (x) = F j (x ) because both the input to node i and the input to node j change between x and x . In this case the STG is not well defined on the extended space X.
We resolve this problem by unfolding the parameter p = (θ, f (θ)) into two parameters p 1 = (φ, g(φ)) and p 2 = (π, h(π)) which correspond to splitting the two thresholds ζ i,i = ζ j,i into their two possible orderings ζ i,i < ζ j,i (Figure 2c) and ζ j,i < ζ i,i (Figure 2d), respectively. The unfolded parameters p k for k = 1, 2 are parameters for the network
which is the transport network with the edges e 1 and e 2 no longer identified as a transport pair and removed from E tr . The unfolded order parameter agrees with the original parameter on each node u = i where we require φ u = π u = θ u . Let * := θ -1 i ({e 1 , e 2 }) be the rank of the transport edge pair {e 1 , e 2 } in TN according to the order parameter θ i . For node i, each new order parameter for the network TN
preserves the ordering given by θ i
and specifies a a particular ordering between edges for e 1 and e 2 :
ζ ii (a) The unfolded logic parameters g(φ) and h(π) agree with the original parameter on each node u = i, i.e. g u = h u = f u . For node i, an additional function is added to reflect the additional threshold: g i := (g q ) q∈{e 1 ,e 2 }∪T r (i)∪(T tr (i)\{{e 1 ,e 2 }}) and h i := (h q ) q∈{e 1 ,e 2 }∪T r (i)∪(T tr (i)\{{e 1 ,e 2 }}) .
These functions satisfy g q = h q = f q for q / ∈ {e 1 , e 2 } and
This implies that the logic parameters for both unfoldings agree g = h.
The asynchronous Boolean dynamics for the parameter (θ k , f k (θ k )) are then defined on the expanded STG X which contains an additional state in the ith direction
If the resulting dynamics have an inconsistent black wall, then the unfolding procedure is repeated until all inconsistent black walls have been removed.
We view the unfolding as a process of embedding of a degenerate systems into a family of nearby non-degenerate systems. The switching where two thresholds coincide is degenerate, while switching systems where all thresholds are distinct is generic. At the same time, we view switching system as a degenerate limit of a family of smooth nonlinearities [31, ?, ?]. The correspondence between Morse graphs of the switching systems and Morse sets of the nearby smooth dynamics has been established in some cases in [15,14,22] which shows that dynamics of any unfolded switching system perturbs into nearby smooth systems.
We illustrate our work on transport regulation on a simplified model of p53-Mdm2 regulation.
4.1. p53-Mdm2 network. Protein p53 is a transcription factor that regulates genes involved in growth arrest, DNA repair, apoptosis and cellular senescence [21,35,40], and thus plays a key role in the control of the cell proliferation. Normally, the level of p53 remains low. It becomes elevated only when cells are stressed or damaged. The importance of p53 led to high level of interest in studying the varying response of p53 to different types as well as severity of stress. In particular, Geva-Zatorsky et.al. [25] observed great variability in oscillatory response of individual cells to DNA damage. The oscillation period decreases with increasing irradiation intensity, with a typical cycling time of about 10 hours at low irradiation doses and about 6 hours at high irradiation doses.
They also observed that in some cells, the oscillations stopped or changed frequency after one or two days. To explain some of the diversity in response, Ouattara, Abou-Jaude and Kaufman [1,36] developed a simple 3-dimensional logical model of p53 function incorporating the function of the ubiquitin ligase Mdm2. Their logical bifurcation analysis suggested that the system is capable of birhythmicity, that is, two coexisting stable periodic cycles of different frequencies, for a certain range of irradiation dose [1]. This has been confirmed by bifurcation analysis of a closely related nonlinear differential equation model with steep Hill type nonlinearities [36]. Their logical analysis suggests that there is a set of domains in the state transition graph that form a figure eight graph with two distinct, but overlapping loops, see Figure 3a.
In a follow-up paper Abou-Jaude, Chaves and Gouze [2] used switching function approximations of Hill nonlinearities (7) and studied the differential equation model
The variables are p53 (P ), nuclear Mdm2 (M n), and cytoplasmic Mdm2 (M c). The functions s + and s -are the aforementioned switching functions, V r is a compartment volume ratio, the k values are scalar constants representing reaction and production rates, and the d values are decay rates.
The motivation for this model is the following (see [1] for a more detailed description). Nuclear Mdm2 down regulates p53 by increasing its degradation rate through ubiquitination [5,17,29] and by blocking its activity [6,33,34]. On the other hand, p53 up regulates the cytoplasmic level of Mdm2 by increasing the transcription of gene MDM2 [4,19]. Finally, p53 inhibits the transport of
Mdm2 from the cytoplasm to the nucleus [28]. These interactions are summarized in the conceptual network model in Figure 3b that we call the influence graph.
Abou-Jaude, Chaves and Gouze [2] used a 2-dimensional switching system as an approximation of system (7) and showed that there is a region of parameter space that admits birhythmicity. Based on experimental data that measures the abundance of M n and P , the projection of the corresponding periodic orbits have the approximate shape exhibited in Figure 3a.
To adapt the model (7) to our notation, we note that functions s ± in (7) (see [2]) are piece-wise constant with values 0 and 1, while our functions σ ± have values L < U that are free parameters. Therefore we set
and with this notation we arrive at a model
In addition, we make two more substantial modifications of this model. First, we note that the effect of M n on P is inactivation by ubiquitination, which represents a bi-molecular reaction where P must be present to be inactivated. We therefore replace the term σ -(M n) in the last equation by the term σ -(M n)P . Second, in agreement with the DSGRN modeling approach described in this paper, we replace the transport term σ -(P )M c in the second and third equations, which is linear in the transported substrate M c, by the transport term σ -(P )σ + (M c).
We arrive at the following set of equations that we will study using DSGRN dM c dt
Comparing to equation ( 5), the volume ratio V r in the this equation is represented by the parameter A in (5).
There is an apparent inconsistency in our model. We replace the term σ -(P )M c that is linear in M c by σ -(P )σ + (M c) in the first two equations in ( 8), but leave a similar term σ -(M n)P unaltered in the third equation. The reason is that the ubiquitination term σ -(M n)P can already be modeled in DSGRN by writing it as a modification of the decay term of P [11]. This approach does not work for the transport term and therefore we make a modification or the transport term that permits DSGRN-type analysis.
We will search the essential DSGRN transport parameter graph for STGs that admit a stable set of states A that are consistent with the two types of experimentally observed periodic orbits in Figure 3a. That is, the set of states A should admit two periodic paths that in the projection to M n, P components form both a shorter and a longer closed path through the STG. We will call such a set of nodes a figure 8 attractor A 8 . The projection of A 8 into M n, P coordinates should resemble Figure 3a, or its mirror image along the vertical axis.
Since our method is coarse, the existence of A 8 does not guarantee that there are parameter sets where both stable periodic orbits coexist in system (8), as was shown in [2] for their model. Such coexistence may exist only in a small parameter region in the mathematical model and biologically may not be as relevant as existence of both types of orbits for parameters that are near each other. Uncertainties inherent in the modeling process produced by uncertain choice of parameters and nonlinearities, combined with ubiquitous noise in biological system and their experimental measurements suggest that the existence of both orbits at nearby parameters may be sufficient to conclude that the model is capable of reproducing the experimental data. In a DSGRN framework, nearby parameters are defined to be real-valued parameter sets (L, U, ζ, γ) for an ODE system that are consistent with the logic and order parameters at the same DSGRN parameter (θ, f (θ).
Network. Following the development in Section 3.1, we transform the conceptual p53-Mdm2 model in Figure 3b to the network structure in Figure 3c via (8). There are two transport edge pairs in Figure 3c, the pair emanating from M c and the pair emanating from P . Notice that in addition there are multi-edges from P to M c, one that controls transcription and one that controls transport. We will see how these structures impact the formation of admissible sentences derived in Section 3.3. The list of states of X and the corresponding input sequences generated at those states under two order parameters θ and φ are summarized in Table 2. These sequences form a sublattice with partial order induced by order 0 < 1 on all components. In Figure 4a is depicted set I M c (θ) of realizable input words in partial order (1) and in Figure 4b is depicted set I M c (φ) of realizable input words in partial order (2). Since the node P affects node M c via two edges e M c,P and e M c,P the order of these edges induced by the order parameter θ or φ restricts the realizable words.
To define the state transition graph, we count the number of thresholds corresponding to the targets of each node
The key to the definition of a map F 0 is evaluation f e (b j (x; θ)) of logic parameter f e at a Boolean sequence b j (x; θ) at each domain x ∈ X. The transition from the domain (0uv) to (1uv) changes two out of three input sequences (b M c , b M n , b P ): the sequence b M c changes from (1ab) to (0ab) and b M n changes from (1b) to (0b). These simultaneous transitions lead to inconsistent black wall for some parameters (θ, f ). To address this issue, we unfold the parameter p = (θ, f ) to parameters p 1 = (φ, g) and p 2 = (π, h). We use p 1 = (φ, g) to denote the order parameter which orders the self-edge e tr M c,M c first, and use p 2 = (π, h) for the order parameter that orders the self-edge second:
(000) (001,10,0) (000) (001,10,0) (001,10,0) (100) (000,10,0) (001,11,0) (100) (000,11,0) (200) (000,11,0) (000,11,0) Table 3. the STG denoted by their signatures in the order (M c, M n, P ). ( 1) Two neighboring states (000), (100) in X that differ in the first component. The unfolded logic parameter g is then a collection of functions
) . The unfolded logic parameter h is the same as g: h = g.
There are no inconsistent black walls for the dynamics on Y with parameters (φ, g) or (π, h) because the change of b M c from (1ab) to (0ab) happens on a separate threshold from any other changes (see Table 3). Each of these parameters will admit a consistent black wall.
4.2.4. Essential parameter graph. We are interested in description of essential parameter graph that supports existence of attractor A 8 . We therefore assume that none of the logic parameters are constant, see section 3.5 and thus each logic parameter has the value zero for the lowest input sequence, and the value one on the highest input sequence in the partial order.
• The logic parameter
• The logic parameter f P = (f M c,P , f tr(P ) ) satisfies
Theorem 4.1. Assume ( 11), ( 12) and ( 13).
(1) The essential parameter graph P G for system (8) is depicted in Figure 5. The front eight nodes have order parameter θ, the back eight nodes have order parameter φ.
(2) Every bold parameter node admits at least one unfolding which contains A 8 attractor in the corresponding STG.
Parameter graph for both P threshold orderings. In the front are parameter nodes with order parameter θ, in the back with order parameter φ.
Before proceeding with the proof of the Theorem, we discuss some properties of the STG implied by assumptions ( 11), ( 12) and ( 13).
Lemma 4.2. The assumptions ( 11), ( 12) and ( 13) imply the existence of the arrows in Figure 6a.
Proof. We first discuss the structure the STG. There are three levels of STG in the P direction: bottom (B) corresponding to states (M c, M n, P ) = (ab0), middle (M) corresponding to states (ab1), and top (T) corresponding to states (ab2) for a, b ∈ {0, 1}. We first note that since the effect of M n on P is mediated by a function σ -(M n), by definition of b P , the value of b P = 1 in any state (M c, M n, P ) = (a0b), a, b ∈ {0, 1}. Therefore the assumption (13) implies that in all states (a0b) the transitions in STG go in positive P direction, i.e. upward. This determines vertical arrows on the left side of Figure 6a. Similar argument shows that at all states (a1b), a, b ∈ {0, 1}, on the right side of STG in Figure 6a, the transitions go downward.
The assumption (12) implies that in states (M c, M n, P ) = (100) and (110) where the input b M n to f P,M n is 11 the arrow goes from left to right; in the states (002) and (012) where the input is 00, the arrow goes right to left.
Finally, assumption (11) implies that at the states (M c, M n, P ) = (002) and (012) where the input b M c to f tr(M c) is 111 the arrows go forward and at (100) and (110) where the input to f tr(M c) f M c,M c = f M n,M c is 000 the arrows go backward. This completes description of STG in Figure 6a.
To prove the Theorem we separate proofs for the two different choices of order parameter.
Proposition 4.3. Assume ( 11), ( 12) and ( 13).
(1) For order parameter θ all values of the logic parameter f tr(M c) listed in Table 4, except A 1 , have an STG that admits A 8 .
(2) For order parameter φ all values of the logic parameter f tr(M c) listed in Table 5, except A 2 , have an STG that admits A 8 .
The labels of the parameter nodes correspond to positions of black wall in the M c component (front-back direction in Figure 6a) in STG.
There are three levels of STG in the P direction: bottom (B) corresponding to nodes (ab0), middle (M) corresponding to nodes (ab1), and top (T) corresponding to nodes (ab2) for a, b ∈ {0, 1}. The label of the parameter nodes denote the position of the black walls in the corresponding STG; for instance T M B parameter node has black walls on all three levels, while E does not have any black walls. We will show below that the parameter nodes A 1 , A 2 cannot admit A 8 .
Remark 4.4. Since our emphasis is on finding parameter nodes that support A 8 in STG, we only consider essential parameter nodes, see assumptions (11), ( 12) and ( 13). Abou-Jaoude et.al. [1] were interested in collection of STGs as a function of the decay rate of M n. This is equivalent to a path in parameter graph P G(M n) that also includes constant logical parameters f P,M n (00) = 1 and f P,M n (11) = 0 which violate assumption (12). It is easy to see that when f P,M n (00) = 1, the corresponding STG has attracting nodes (110) or (010) and if f P,M n (11) = 0 then the corresponding STG has the attracting states (102) or (002). These observations agree with those of [1] as these attracting states represent reset state after DNA damage repair. Perturbation of parameters to those that support oscillations represent response of the system to DNA damage and pulse generation. 000 001 010 011 110 111 0 Note that this means that the edge e M n,P is input in-essential since the value of the logic parameter f P,M n only depends on the first input component that correspond to σ + (M c). The value of f P,M n does not depend on the second component of the input that corresponds to σ -(P ).
Proof. Referring to Table 2 we note that for states x = (M c, M n, P ) = ( 102 Similarly, note that for states x = (000) or x = (010), the value of b M n (x) = 10 for any choice of order parameter. If f P,M n (10) = 1 then the STG would have an arrow from (000) to (010) which, together with assumption (12) would imply the A 8 attractor does not exist. Therefore we must have f P,M n (10) = 0.
Corollary 4.6. Any STG that satisfies assumptions (11), ( 12) and ( 13) and contains A 8 has to have directions depicted in Figure 6b.
The dotted directions will depend on a choice of a parameter node in Figure 5 and will be discussed next.
We first show the tight correspondence between STGs under order parameters θ and φ, Let α denotes a word substitution ( 14) α(101) = 011, α(100) = 010.
Note that with exception of rows corresponding to A 1 , E 1 and A 2 , E 2 , in all other rows denoted by the same symbol (i.e. B, M, T, T M, M B, T M B) the entries in Table 4 and Table 5 are the same when the column labels are mapped by α.
The same applies to Table 2; the values of all three b M c , b M n , b P between second and third column are identical when we apply substitution α. Therefore the identically labeled rows of tables Table 4 and Table 5 in Proposition 4.3 encode the same multi-valued maps, the state transition graphs will be identical for identically labeled rows. Therefore it is sufficient to describe STGs only for order parameter θ.
We discuss all possible parameter nodes of the factor parameter graph corresponding to M c, each of which is given by a row in Table 4. To construct an STG for a given row, we start with the partial STG in Figure 6b, because each arrow there is necessary for A 8 by Corollary 4.6. The only arrows missing from Figure 6b are in the M c direction. To get these arrows, we use Table 2 and locate an STG node in the left column and, since we are working with order parameter θ, identify the value of b M c in the middle column. Then the value of b M c at the top of the column in Table 4 is an input to f tr(M c) whose output is listed in the appropriate row of Table 4. If this value is 0 the arrow points from front to back, that is, in the negative M c direction. If this value is 1 the arrow points from back to front in the positive M c direction. If f tr(M c) (b M c (x)) = 0 at a state x = (1ab) and f tr(M c) (b M c (x)) = 1 at the adjacent state x = (0ab), then the STG has a black wall. Since STG arrows in Figure 6b point to the right on the front states (1ab) but point to the left on the back states (0ab), these black walls are inconsistent (see Section 3.7 and Figure 2). We discuss the unfolding of the parameters B 1 , M B 1 , T M B 1 , M 1 , and T 1 which produce an inconsistent black wall in the following section.
We now discuss remaining logic parameters f tr(M c) which are distinct from B 1 , T 1 , M 1 , T B 1 , M B 1 and T M B 1 and still consistent with the monotonicity condition (6). These correspond to A 1 :
We observe that the logic parameter A 1 is not compatible with Lemma 4.5. To see this note that
Since the set of realizable words for M c and M n due to common transport term are dual with each other, and these correspond to second and third input, this implies that f P,M n (10) > f P,M n (01), which contradicts Lemma 4.5. Therefore logic parameter A 1 is not compatible with STG with A 8 .
The STG for logic parameter E 1 is depicted in Figure 8a; there are no black walls in STG. The logic parameter of E 1 is in the second row of Table 4.
The argument for order parameters A 2 and E 2 are analogous; in these case, however the word substitution α maps logic parameter E 1 to A 2 , which is not compatible with Lemma 4.5. The substitution α maps A 1 to E 2 which is in row 6 of Table 4 and its STG is in Figure 8b.
To complete the description of the parameter graph depicted in Figure 5, note that the nodes with subscript 1, that are connected by edges are related by logical adjacency, since the corresponding rows in Table 4 differ by exactly one value. The nodes connected between back and front walls which differ by the subscript, are connected by order adjacency. so that superscript R represents the back ("rear") unfoldings and superscript F represents the front unfoldings.
We start with extending partial STG in Figure 6a.B by adding arrows in the Mc direction, see Figure 7. These arrows depend on the parameter chosen. At each parameter in Figure 7, there is a black wall between at least two states of the partial STG. These black walls are denoted by a black dot. Note that the name of the parameter indicates the position of these black walls in the P direction: B, M and T stands for bottom, middle and top black wall. The transitions along the inconsistent black walls in the Mn direction (i.e. left to right) are not well defined as they can either agree with the arrows on the back wall (right to left) or on the front wall (left to right). In each partial STG the red dashed arrows are hypothetical transitions which, if realized, would complete A 8 attractors together with the other red arrows. We now determine if these hypothetical arrows give rise to actual transitions in an STG in at least one unfolded parameter.
In both unfoldings the black wall, which is a subset of ζ M c,M c , becomes consistent, since the M n directions on both sides of ζ M c,M c will now agree. Therefore under either unfolding there is a well defined state transition in the M n direction (either to the left or to the right in the STG) along the black wall. The STG with either unfolding gains a new state that represents the new set of domains bounded between the two thresholds. We denote these new states as empty circles in Figure 9, which shows the STGs for the parameters and unfoldings which admit A 8 . As an example, compare the partial STG at parameter B 1 in Figure 7a with STG for the back unfolding of B 1 in Figure 9a. Note that the dashed hypothetical transition in Figure 7a is replaced by a regular transition in Figure 9a and that these red transitions form A 8 . To continue the argument, observe that in order to turn hypothetical dashed arrows that go from right to left in partial STGs of parameters B, M B and M (Figures 7a, 7b and7d), the M n direction along the black wall must coincide with the M n direction along the back side. This implies that
which is realized by back unfolding, see Figures 9a, 9b and9c.) The same argument shows that the front unfolding turns hypothetical dashed arrows from left to right in partial STGs of nodes T, T M and M (Figures 7f, 7e and7d) into actual transitions in STG, see Figures 9f, 9e and9d.) Note that either unfolding works for parameter node M since in both cases there is an alternative path between the front and back walls of STG. On the other hand, neither unfolding works for T M B because front unfolding prevents transition from states (112), ( 111) and (110) towards the back wall, and back unfolding prevents transition from (002), ( 001) and (0) towards the front wall. To summarize,
• Back unfolding of STG at logic parameters B and M B admits A 8 , Figures 9a and9b.
• Front unfolding of STG at logic parameters T and T M admits A 8 , Figures 9e and9f; • Both front and back unfolding of STG at logic parameters M admit A 8 , Figures 9d and9c, respectively. We will use notation M F 1 and M R 1 to denote front and back ("rear") unfolding of parameter node M 1 , respectively.
• Neither unfolding of STG at logic parameters T M B admits A 8 .
As discussed in section 3.7, switching system with coinciding thresholds is result of taking two limits; one is the limit in which thresholds are converging together, and one is the limit of smooth nonlinearties that converge to a switching nonlinearity. Therefore, as suggested by [15,14,22] if there is a correspondence between the A 8 in the switching system and a system with sufficiently smooth nonlinearities, then even if only one unfolding admits A 8 , then there is an open set of unfolding parameters (distance between thresholds, steepness) around (0,0) that will exhibit A 8 .
To illustrate effectiveness of DSGRN to quickly scan the parameter space and find interesting dynamics like attractor A 8 , we show that the presence of A 8 succesfully predicts existence of both large and small (defined below) stable periodic orbits for an associated set of ODEs with high Hill coefficient. We numerically simulate the following set of ODEs: Note that all the parameters of this system, except the Hill exponent n, are in direct correspondence with parameters of (8). Therefore we fix n = 90 and then we can associate each parameter set of (15) to a particular parameter node of DSGRN model (8). Since the parameter node is given by a set of inequalities, we sample each of them and then compare results of ODE with predictions of DSGRN offered by presence of A 8 . We employ Mathematica's FindInstance function to generate one parameter set within each of the unfolded nodes
and T F 1 , where the superscript indicates the type of the unfolded parameter node. Using this parameter set to determine the L and U values, we then uniformly sample 100 values for the thresholds ζ within the range determined by the parameter node and these L and U values. The resulting samples are not expected to be representative of the parameter node.
6.1. Data analysis. For each parameter set we generate 10 initial conditions in the phase space, chosen in the DSGRN domain which corresponds to the node in the STG with two different out edges. This is done in order to maximize the probability that some of the solutions will follow one path and others a different path in A 8 (e.g. node 111 in Figure 9c). In each version of A 8 attractor we identify small and large loop, see Figure 7. Since every loop visits both small and large values of M c, and M n, the distinction is in the P variable. The large loop crosses all three levels of P while small visits only two of them. For example, see Figure 9d. The small loop visits only the lower two levels of P values and the large loop visits all three levels. While there are three distinct paths that correspond to large loop, we do not distinguish between them for this classification.
For a particular parameter sample, we numerically integrate the ODE. For all sampled parameters and all initial conditions, the numerical solution converges to a periodic orbit. We classify each limit periodic orbit as either a small or a large periodic orbit, i.e. one that lies within the small vs. large loop of A 8 . Since the A 8 is an attracting region for system (8), but not necessarily for (15), we classify an orbit as large if the orbit crosses both P thresholds and the orbit is small otherwise. See Appendix A for representative examples of large and small orbits found by numerically integrating (15). There are two levels at which birhythmicity can be established. On one hand, two stable periodic orbits can coexist at the same parameter set for (15). This is the behavior observed for a very narrow range of parameters in [2]. This would manifest itself in our simulation as a presence of parameters for which (15) has some initial conditions converge to a large periodic orbit and some initial conditions converge to a small periodic orbit.
A coarser, yet perhaps biologically equally relevant behavior, is when at one parameter value (15) exhibits a small periodic orbit for all initial conditions, while for another nearby parameter value (15) exhibits a large periodic orbit for all initial conditions. For the purpose of this second approach we consider parameter values near each other if they are sampled from the same DSGRN parameter node.
, we have identified some parameter samples which yield small orbits and other parameter samples which yield long orbits. This verifies that the second, coarser version of birhythmicity exists in the ODE model (15). The proportion of samples resulting in small versus large orbits appears to be parameter node-dependent (see Table 6). However, we do not expect that the relative frequencies of small and large orbits reported here to reflect the true frequencies across the whole parameter node because our samples are not representative of the whole parameter node. Rather, we interpret these results to mean that most parameter sets within a parameter node will exhibit stable periodic orbits and both small and large orbits are possible.
Having found both large and small orbit parameter sets in each parameter node, we employed the following strategy to search for a parameter set for which both a large and small periodic orbit coexist. Given a parameter node, we selected a pair of parameter sets within the parameter nodeone parameter set corresponding to a large orbit and one parameter set corresponding to a small Parameter node Small Large Table 6. Summary of behavior for 100 parameter samples within each parameter node. Counts of parameter samples which yield small or large orbits for all initial conditions tested.
orbit. We then drew a line in parameter space connecting these two parameter sets. Along this line we expect either that the large orbit smoothly transitions into the small orbit or that there is an interval in which two orbits coexist. To examine which transition occurs, we sampled equidistant parameters along this line. We repeated this process across all parameter nodes and all large-small parameter pairs within each node. In all cases, we found that the periodic orbits transitioned smoothly between large and small (data not shown).
The language of networks is increasingly used across scientific disciplines. Undirected networks often reflect correlations or similarity between entities represented by the nodes, while the directed networks aim to represent causality, or relationship in time. Directed networks have been used in systems biology to describe effects between genes, proteins and small molecules like metabolites. The description of the dynamics of the molecular interactions, represented by directed networks, has been a focus of mathematical biology for a long time. These models can use differential equations, Boolean formalisms or stochastic descriptions, but in all of these modeling approaches it is important to realize that the edges of the network may represent different biological phenomena that require different models.
In this work we investigate how to incorporate transport regulation into the combinatorial modeling approach DSGRN [10,24] that was originally designed to model transcriptional regulation and only recently extended to include post-transcriptional modifications [11]. DSGRN is able to describe all potential dynamics of the network by imposing the strong assumption on network interactions that they are switch-like. This has the advantage of achieving a large, but finite and thus directly enumerable description of the dynamics, that, at the same time, provably describes the dynamics of collection of ODEs with steep nonlinearities [22]. The mathematical structures that enable such a finite enumeration of dynamics by DSGRN are not well suited for transport regulation modeling, where a compound regulates rate of transfer between two compartments. The results of this paper show that by generalizing the DSGRN approach, one can include transport regulation in the class of cellular processes that can be modeled by DSGRN.
Comparison of models of cell regulation to experiments is often difficult because of presence of measurement noise, ill-defined cellular conditions and necessary simplicity of the models. The difficulty with parameterization of traditional ODE models, as well as uncertainty in the process of selection of types of nonlinearities used, makes traditional comparison using differences between a particular solution (i.e. a trajectory of the model and often only few experimental data points) problematic. An alternative to this trajectory based approach is comparing models and experiments on a more qualitative way. Is the model capable of producing bistability? Can this model support a stable oscillation? These comparisons are then used to reject particular models or networks rather than selecting a unique model capable of matching the data. This process reduces hypothesis space which can be very useful to the experimentalist and avoids making conclusions based on a narrow fit of a particular model, at a particular parameter and for a particular initial condition.
We illustrate our extension of DSGRN to transport modeling by applying our method to a p53-Mdm2 network model. The experimental observations show the existence of two types of stable oscillations. This had been previously interpreted as co-existence of two stable periodic orbits [1,36] and indeed, such co-existence had been found in two different models for a particular and narrow range of parameters [2]. Following our philosophy of seeking a qualitative agreement with the experiments, we look for the existence of such a pair of orbits in the model at two potentially different, but close, parameter values. We reason that during an experiment, especially in human cells, such shifts in parameters occur all the time and are not experimentally detectable. With this lower standard for birhythmicity, we found that birhythmicity is common across parameter space and verified this DSGRN prediction by simulation of the associated smooth ODE model.
Parameter Fig.