Control of Intracellular Molecular Networks Using Algebraic Methods

Many problems in biology and medicine have a control component. Often, the goal might be to modify intracellular networks, such as gene regulatory networks or signaling networks, in order for cells to achieve a certain phenotype, such as happens in cancer. If the network is represented by a mathematical model for which mathematical control approaches are available, such as systems of ordinary differential equations, then this problem might be solved systematically. Such approaches are available for some other model types, such as Boolean networks, where structure-based approaches have been developed, as well as stable motif techniques. However, increasingly many published discrete models are mixed-state or multistate, that is, some or all variables have more than two states, and thus the development of control strategies for multistate networks is needed. This paper presents a control approach broadly applicable to general multistate models based on encoding them as polynomial dynamical systems over a finite algebraic state set, and using computational algebra for finding appropriate intervention strategies. To demonstrate the feasibility and applicability of this method, we apply it to a recently developed multistate intracellular model of E2F-mediated bladder cancerous growth, and to a model linking intracellular iron metabolism and oncogenic pathways. The control strategies identified for these published models are novel in some cases and represent new hypotheses, or are supported by the literature in others as potential drug targets. Our Macaulay2 scripts to find control strategies are publicly available through GitHub at https://github.com/luissv7/multistatepdscontrol.

However, increasingly many published discrete models are mixed-state or multistate, that is, some or all variables have more than two states, and thus the development of control strategies for multistate networks is needed. This paper presents a control approach broadly applicable to general multistate models based on encoding them as polynomial dynamical systems over a finite algebraic state set, and using computational algebra for finding appropriate intervention strategies. To demonstrate the feasibility and applicability of this method, we apply it to a recently developed multistate intracellular model of E2F-mediated bladder cancerous growth, and to a model linking intracellular iron metabolism and oncogenic pathways. The control strategies identified for these published models are novel in some cases and represent new hypotheses, or are supported by the literature in others as potential drug targets.
Keywords Intracellular network · Discrete dynamical system · Control · Polynomial dynamical systems 1 Introduction Modification or differential regulation of intracellular networks, either at the gene, protein, or metabolite level, can result in an altered cellular phenotype. Being able to perform targeted modifications for this purpose may be desirable for several reasons, whether to alter bacterial metabolism for industrial production [38] or to mitigate properties of a tumor cell [16,17,59]. Systematic approaches to the identification of such targeted modifications are therefore of considerable importance. Generally, this is accomplished through the use of mathematical models as discovery tools. In addition to systems of differential equations, an increasingly common modeling framework are time-and state-discrete models, such as Boolean networks and their various generalizations. These provide only semi-quantitative information but are more easily constructed, since they do not require quantitative kinetic information, and they can sometimes be more intuitive for the experimentalist. A drawback of discrete models is that their underlying mathematical theory, in particular control theory, is not yet well-developed, and the present paper makes a contribution to this body of work. Given a mathematical model of an intracellular regulatory network, one commonly associates the possible phenotypes of the cell with the attractors of the model, an idea that can be traced back to Waddington [50,42] and Kauffman [23,21]. For example, the steady states in [37], discussed in more detail below, correspond to proliferative, apoptotic, or growth-arrest phenotypes of a cancer cell. In [3], the steady states of the model correspond to the observed altered iron metabolism phenotypes in a breast epithelial cell with and without a certain RAS mutation. We present here a method to systematically identify modifications to a model that can change the attractor landscape in prescribed ways. Such modifications consist of, e.g., deactivating a node or modifying the effect of an edge in the model's wiring diagram, a graph-theoretic representation of the functional dependencies of the different model variables. In the case of gene regulatory or signaling networks, this could be accom-plished, for instance, through a compound that blocks the protein corresponding to a particular gene.
We focus on the mathematical modeling framework of multistate discrete networks. These can be defined as dynamical systems that are discrete in time as well as in variable states. More formally, consider a collection x 1 , . . . , x n of variables, each of which can take on values in a finite set X 1 , . . . , X n . Let X = X 1 × · · · × X n be the Cartesian product. A discrete dynamical system in the variables x 1 , . . . , x n is a function where each coordinate function f i : X → X i represents how the future value of x i depends on the present values of all the variables. If X i = {0, 1} for all i, then each f i is a Boolean rule and F is a Boolean network.
Discrete networks defined in this way can be represented in the richer mathematical framework of polynomial dynamical systems (PDSs) [48,19], as can models in other common frameworks, such as Boolean networks [1], logical regulatory graphs [2], or multistate networks [37,3,11,45]. For instance, the mathematical tools associated with PDSs allow for the computation of all steady states and cycles up to a certain length of a system as the solutions to a system of polynomial equations (in a suitably chosen finite field) without explicit simulation of the entire state space. Recently, we have used these mathematical tools to construct a method that rigorously computes modifications for the control of Boolean networks that can help avoid regions of the state space or create new steady states or cycles [32]. The method was applied to a mathematical model of cellular response to DNA damage [4] and a model of large granular lymphocyte apoptosis escape [60]. However, an increasing number of discrete mathematical models in this context are not Boolean, e.g., [37,3,11,45], so that a more general method is desirable, and such a generalization is the focus of this paper.
We note that there are several published control methods that do not rely on the PDS representation of discrete models, such as Stable Motifs [55], Feedback Vertex Sets [57], Minimal Hitting Sets [49,25], and several others [36,26,35,14,56]. Our method provides a flexible control framework that, for instance, allows for the identification of controllers for creating new (desired) steady states, a feature that other methods do not allow.
It also extends our method for Boolean networks [32,31] to multistate networks, and thus broadening the scope of use of the PDS representation of discrete models.
As our method uses polynomial algebra over a finite field, all network nodes need to take values in a common finite field, in particular, all nodes need to have the same number of possible values. In many published models, however, different nodes take on different numbers of states, and this number generally does not allow the imposition of a finite field structure (for which the number is required to be a power of a prime number); see, e.g., [37,56]. As part of the algorithm in this manuscript, we present a method to convert models with a general number of mixed discrete states into a model that satisfies the computational algebra requirements, without changing the model's steady states, and which is not equivalent to the well-known reduction to a Boolean network that adds new nodes to the network, as done in [56].
We demonstrate the power and versatility of this method by applying it to two recently published multistate network models. One is a model of bladder cancer response to different stimuli, including DNA-damage, EGFR, FGFR3, and growth inhibitors [37]. The method can find combinatorial interventions that block proliferative steady states. The second is an intracellular iron network model in breast epithelial cells presented in [3]. We identify interventions to recover basal expression of the iron export protein of a malignant breast epithelial cell with RAS over-expression. These can be viewed as predictions to be validated experimentally.

Methods
As the first step of the method, a multistate discrete network is represented in an algebraic framework. In the process, we provide a general procedure to extend any multistate discrete network to a polynomial dynamical system. We refer the reader to the following books [7,27] for the basics of finite fields and the basics of computational algebra.

Discrete dynamical systems
In this manuscript, we consider discrete variables x 1 , . . . , x n , each taking on values in a finite set X 1 , . . . , X n .
For the purpose of exploiting the algebraic properties of discrete functions, it is assumed that the variables x 1 , . . . , x n take on values on a finite field F. Then, using the fact that any discrete function f i : F n → F can be represented as a polynomial on {x 1 , . . . , x n }(see e.g. [22,27]), that is f i ∈ F[x 1 , . . . , x n ], the discrete network can be represented as The discrete network can thus be represented as a polynomial dynamical system (PDS) [48,19]. We describe how to convert a mixed-state model into a PDS in the Appendix, Section 7.1. We also provide a concrete example of how to carry out this procedure when the update rules are written as Boolean expressions of conditions in the Appendix, Section 7.2.
Given a discrete network F = ( f 1 , . . . , f n ), we can define its wiring diagram W to be the directed graph with n nodes x 1 , . . . , x n associated to F, such that there is a directed edge in W from x j to x i if x j appears in f i , and there exist a 1 = (a 1 , · · · , a j−1 , b j , a j+1 , · · · , a n ), a 2 = (a 1 , · · · , a j−1 , c j , a j+1 , · · · , a n ) ∈ F n such that In other words, the value f i takes on depends on the values of x j . The dynamics of discrete networks are given by the difference equation x(t +1) = F(x(t)); that is, the dynamics are generated by iteration of F. More precisely, the dynamics of F are given by the state space graph S, defined as the graph with vertices in F n which has an edge from x ∈ F n to y ∈ F n if and only if y = F(x). In this context, the problem of finding the states x ∈ F n where the system will get stabilized is of particular importance. These special points of the state space are called attractors of a discrete network and these attractors may include steady states (fixed points), where F(x) = x, and cycles, where F r (x) = x for some integer r > 1. Attractors in network modeling might represent a differentiated cell type [24] or a cellular state such as apoptosis, proliferation, or cell senescence [20,41]. Identifying the attractors of a system is an important step towards system control. For example, a steady state might represent a cellular phenotype characterized by low expression of an iron exporter [3], a particularly deleterious phenotype in several cancers [46].

Control actions
In this manuscript, we focus on controlling the dynamics of a multistate network by avoiding undesirable steady states, creating new steady states, or avoiding regions in the state space, accomplished by modifications to the wiring diagram of F. This is an extension of our previous control methods applicable to Boolean networks [32].
We consider two types of control actions: 1. deletion (or constant expression) of edges and 2. deletion (or constant expression) of nodes. An edge deletion represents the experimental intervention that prevents a regulation from happening. These actions can be achieved by the use of therapeutic drugs that target a specific gene interaction [4]. Constant expressions could also help to drive the system into a more desirable state [39].

Encoding control actions in multistate networks.
In the Boolean setting, the deletion of an edge in the wiring diagram was implemented by setting an input to zero so that the interaction of that input (represented by an edge) was being silenced. For the multistate case, the silencing of the interaction will be applied whenever the control variable is within a range of values of the possible discrete values. For expository purposes, we use a simple function taking the value 1 on the singleton set {0} and 0 on the set F − {0}.
Definition 1 (Edge control for multistate networks) Consider the edge x i → x j in the wiring diagram W and let q =| F |.
For u ∈ F, the control of the edge x i → x j consists of manipulating the input variable x i for f j in the following way: For each value of u ∈ F we have the following control settings: . This is the case when the control is active and the action represents the removal of the edge That is, the control is not active.
Similarly, for node deletion in the multistate setting we have the following definition.
Definition 2 (Node control for multistate networks) Consider the node x i in the wiring diagram W and let q =| F |. The function encodes the control of the node x i because for each possible value of u ∈ F one has the following control settings: This action sets the function of x i to zero. For instance, this can represent the knock-out of gene x i or blocking the synthesis of a protein . That is, the control is not active.
We note that in Equations 1-2 we only consider edge and node deletions. In general, one could consider setting edges and nodes to a constant value within F.

Generating new steady states.
Suppose that y 0 = (y 01 , . . . , y 0n ) ∈ F n is a desirable cell state (for instance, it could represent the state of cell senescence in a cancer model) but is not a steady state, i.e., F(y 0 ) = y 0 . The problem, then, is to choose a control u ∈ F such that F (y 0 , u) = y 0 . We now show how this can be achieved in our framework.
After encoding the multistate network with control as a PDS we consider the system of polynomial equations in the u parameters:

Destroying existing steady states or, in general, blocking transitions
Suppose that x 0 ∈ F n is an undesirable steady state of F(x), that is, F(x 0 ) = x 0 (for instance, it could represent a tumor proliferative cell state that needs to be avoided). The goal here is to find a set of controls such that F (x 0 , u) = x 0 . More generally, one may want to avoid a transition between two states x 0 and z 0 . That is, we want to find controls such that F (x 0 , u) = z 0 . To solve this problem consider the following equation,

Blocking regions in the state space
We now consider the case where we want the dynamics to avoid certain regions of the state space. For example, if a particular value of a variable, x k = a ∈ F, activates an undesirable pathway, or is the signature of an abnormal cell state, then we want all steady states of the system to satisfy x k = a. In this case, we consider the system of equations Note that in contrast to previous sections, we are now using variables for x instead of specific values. Since the steady states with x k = a are to be avoided, we want to find controls u for which Equation 5 has no solution.

Identifying control targets
In each case of the tasks above we obtain a system of equations (or a single equation) that we need to solve to find the appropriate controls. This can be done using computational algebra tools. For instance, we can compute the Gröbner basis of the ideal associated with Equation 3, see [6], The computation of a Gröbner basis allows us to read off all controls as the solutions to the system of equations. Furthermore, the algebraic approach can detect combinatorial control actions such as control by the synergistic combination of more than one action; see Section 3 for examples.

Results
We first apply the control methods to the mathematical intracellular bladder cancer model in [37]. Then we will apply our methods to the intracellular iron network model in breast epithelial cells in [3]. As a sample control problem, we first show how we can find control strategies for blocking proliferative steady states by blocking interactions and nodes. Then we will identify interventions to recover a desirable fixed point in a malignant breast epithelial cell with RAS over-expression.

Description of E2F-mediated growth model
In [37], the authors present a generic network centered around how a cell in response to stimuli such as the growth factors FGF3, EGF3, and nodes representing growth inhibition and DNA damage ends up in different states such as a proliferative state, apoptotic state, or growth arrest. The model includes several of the cyclins which regulate cell cycle progression by interactions with cyclin-dependent kinases, as well as  [37]. This figure was generated using GINsim [33]. The red arrows represent inhibitory interactions and the blue arrows represent interactions with an activation/positive effect.
the E2F family of transcription factors, which are released from pRB inhibition in the G1/S cell cycle state and control the transcription of several factors relevant to complete the cell cycle. The model is a multistate model with a total of 30 nodes, where 25 nodes are binary and 5 nodes are ternary. Four nodes serve as input, namely FGF3 , EGF3, Growth Inhibition and DNA damage. Three nodes serve as output, namely Proliferation, Apoptosis, and Growth Arrest; see Figure 1 for the wiring diagram. This mathematical model can then be used to see which inputs lead to a cancerous (proliferative) phenotype, or to generate hypotheses on knockdowns, and/or overexpression of proteins that will evade the proliferative steady state. Due to the size and connectivity of the model, it quickly becomes apparent that predicting emergent behaviors from the blocking of interactions and nodes in the model is not easily done.
3.1.1 Converting the model into a polynomial dynamical system over F 3 In order to use the PDS framework, we first convert the model in [37] into a polynomial dynamical system over F 3 . To do so, we expand the state space to F 3 but in such a way that the steady states do not change. See the appendix for details. The correspondence between variables in the polynomial dynamical system and nodes in Figure 1 is given in Table 1. We used a Docker image of Macaulay2 v: 1.14 [15] and a script available in the Github repository for the conversion.

Avoiding proliferative steady states
In the original model (see Fig. 1), there are four steady states with proliferation equal to 1 (x 24 = 1, Table 1), that is, steady states in which the cells proliferate. It is interesting to note that these steady states have TP53, PTEN, and p21CIP activity at 0, which are all well-known tumor supressors. As mentioned before, steady states correspond to cellular phenotypes. A proliferative steady state thus potentially corresponds to a cell undergoing uncontrolled proliferation. A natural biological question thus arises. How do we avoid such steady states? Our control method predicts several possible control strategies, of which six result in non-proliferative steady states (Table 2). Additional control methods that destroy the original steady states can be combined to block all possible proliferative steady states.   [37] where Proliferation=1.
It is often difficult to generate drugs that target specific proteins, as is often the case of 'undruggable' targets [8]. Moreover, the knockdown of a node could be attained by knocking out a particular gene, but proteins have multiple indispensable physiological actions, and thus a gene knockout might be lethal. It is thus of interest to also be able to affect interactions between two products, which in our context, can be done by targeting edges in the wiring diagram.
We consider the blocking of edges if and only if source and target nodes are both in the set NODES. We also allow node control for nodes x 1 to x 23 . We encode the control of edges and nodes as in Section 2.2.1, and find all solutions to the system encoded by the ideal I = {F − s i } 4 i=1 . It is unfeasible to solve these equations by hand, so we use computer algebra. We compute generators for the Gröbner basis of I, and look at the generators comprised of us alone. The computation of the generators of the Gröbner basis took an average of 3.3653s to compute, with a standard deviation of 0.07636615. We get the following control variables {u 16 , u 14 , u 13 , u 12 , u 10 , u 9 , u 8 , u 3 , u 2 , u 3,10 , u 3,8 , u 2,3 , u 9,16 , u 10,9 u 3,9 , u 11 u 15 }.
In order to avoid the steady states, we can choose any of these generators to not have value 0. The u i s correspond to nodes that are candidates for a knock-down, and the u i, j s correspond to interactions from node x i to x j to block. The product of two control variables indicates that both controls should be applied simultaneously. For example, u 10,9 u 3,9 indicates that the interactions from E2F3 and RAS to E2F1 should be simultaneously inhibited. We present all generated controls in Table 2   It is worth exploring mixed edge-node controls to attempt to destroy all proliferative steady states. For example, Control # 17 in Table 2 is derived by mixing Control # 12 with Control # 5.
Importantly, we see that instead of targeting RAS (Control #3), which has proven to be a formidable challenge and considered a holy grail of cancer therapeutics [43,44], deleting the interaction between RAS and SPRY (Control #13) or FGF3 and RAS (Control #14) has the same effect (in terms of number of proliferative steady states and total steady states) as targeting RAS. Importantly, as RAS is the only node activating SPRY in this model, inhibiting the interaction between RAS and SPRY (Control #13) is equivalent to silencing SPRY, which has been shown to be beneficial in a xenograft model of rhabdomyosarcoma tumors [40]. Furthermore, it has been suggested that the FGF3 mediated RAS activation (Control #14) can lead to Vemurafenib resistance in melanoma cells [51]. Interestingly, inhibiting the CDC25 (Control # 9) family has been suggested as a potential therapeutic for triple negative breast cancer [28], and we see that knocking down CDC25A results in no proliferative steady states.
We thus see that we find some control strategies that destroy all proliferative steady states, some of which have shown promising results in different types of cancers. We also see that some control strategies are not completely effective (for example, knocking down E2F1 alone (Control #5) resulted in four new proliferative steady states, Table 2). However, the control candidates presented here quickly narrow down our list of possible targets. Furthermore, by using function composition, we can also set up control strategies to avoid cycles of a given length.

Algebraic methods applied to the Booleanized network model
The original model presented in [37] is naturally presented in a Boolean network model. It is natural to wonder how the algebraic methods of control compare if we apply them to the representation over F 3 or to the original representation over F 2 . We applied the same methods of control to the network as explained in Section 3.2 (code to do this is in the Github repository). The computation of controls applied 10 times in a Dockerized image of Macaulay2 to the Booleanized model took an average of 6.3528 s with a standard deviation of 0.1718745 (including the running and removal of the Docker container), which is almost twice the time we had for the multistate case (AV=3.37 and SD=0.08). Notice that for the Booleanized network we add an extra node for each multistate node to represent the different possible levels. For example, E2F1 now becomes E2F1:1 and E2F2:2 to represent the possible values. In general, the number of variables to be added could rapidly grow and add computational complexity to the Groebner basis computation. Moreover, without Booleanization, the control nodes and edges are directly related to the wiring diagram of the biological system, making the results easier to interpret and more readily actionable.

Comparison to the Feedback Vertex Set (FVS) control method
Methods of control based on the structure of the network have been developed by several groups. For example, Zañudo et al. [58] proposed to control a set of nodes intersecting all feedback loops in the network (Feedback Vertex Set) plus the source nodes of the network to attain controllability of the network (other groups have suggested using the feedback vertex set as a control target [30,12]). We used the code provided in https://github.com/yanggangthu/FVS_python to approximate a minimal Feedback Vertex Set presented in Figure 1. This yielded the set of nodes {TP53, E2F1, EGFR, RB1, GRB2, CyclinE1} which is a much larger set of nodes to control. It should be remarked however, that the control goals presented in [58] are more general than our narrow control goal of destroying the existing proliferative steady states.
The method presented in this manuscript allows for more targeted and specific control strategies since we are taking the dynamics of the network into consideration. For example, we encoded our problem as destroying the existing proliferative steady states in which case, we observed control sets of size one. Our method can also be implemented as avoiding solutions to the system {F − X, Proliferation − 1} to find controls such that no proliferative state exists.
Some methods of control using the dynamics of networks include methods such as stable motifs for guiding the network towards desired attractors or away from undesired attractors [54,55], using the concepts of Boolean canalization, [5,31], and using the concept of the logical domain of influence [52]. To the best of our knowledge, none of these methods have been implemented for general multistate systems. Interestingly, the concept of stable motifs has recently been generalized for multistate networks [13].

Recovering the iron export capability of breast epithelial cells
Iron is an essential metabolite for eukaryotic cells. Iron is necessary for heme biosynthesis, iron-sulfur cluster generation, and acts as a co-factor in several cellular processes such as DNA replication. It is well accepted that iron metabolism is deregulated in several cancers [46,29]. In particular, low expression of ferroportin, the only known non-heme associated iron exporter in mammalian cells, is associated with poor prognosis in breast cancer [34]. In [3], the authors present a ternary (each variable has three states) mathematical model of how iron metabolism interacts with oncogenic pathways (See Figure 2 and Table 3).
The model presented in [3] predicts that over-activation of RAS leads to low expression of the iron exporter ferroportin. As previously remarked, low expression of ferroportin correlates with poor prognosis. Hence, regaining basal ferroportin levels might prove beneficial to reduce cell proliferation and, thereby, tumor growth, as observed in mice [34].
We thus encode our control problem as avoiding steady states where ferroportin expression is low. To do this, we allow downregulation of the nodes that are not part of the oncogenic pathway (RAS, SOS, ERK, c-MYC, GAPs, EGFR). In other words, we allow downregulation of the following nodes.
We encode the control problem as avoiding solutions to the system encoded by the ideal (F − X, x 3 ) where x 16 = 2 (RAS is overexpressed, see Table 3). We note that 10 computation of the generators of the ideal I = (F − X, x 3 ) took an average of 10.5875s with a standard deviation of 0.2904297, including the overhead of starting the Docker Macaulay2 image.
One of the generators is u 6 * u 15 , that is, simultaneous knock-down of IRP2 and interleukin-6. After setting u 6 , u 15 to one, that is knocking down IRP2 and interleukin-6, and recomputing the steady states, we get exactly one steady state with ferroportin =1, as desired. That is, our model predicts that knocking down IRP2 and interleukin-6 induction of hepcidin will restore the ability of the cell to export iron. This leads to a hypothesis of whether recovering the ferroportin expression of cancerous cells by knocking down interleukin-6 and IRP2 could lead to cell cycle arrest, similar to ironchelators [53]. In fact, ferroportin overexpression in prostate cancer cells has shown some interesting results such as induction of autophagy and p21 overexpression [9].   It should be remarked that the original model of [3] contains a continuity condition, e.g. a state cannot change by more than one unit at each time step. In the appendix, Section 7.3, we show that we can remove this condition if the only thing we are interested in is steady states. In this way, it is easier to interpret controls straight from the interaction network.

Discussion
Encoding discrete dynamical systems as polynomial dynamical systems (PDSs) offers a rich toolbox for the analysis of such models. For example, the computation of steady states and cycles can be encoded as a computation of all the solutions to a system of polynomial equations over a finite field [48], and does not require simulations. Moreover, any discrete dynamical system as defined in this manuscript can be encoded as a PDS, and thus the PDS framework offers an encompassing and general framework for encoding discrete models of biological systems. We previously showed how tools of computational algebra can be applied to find control strategies in Boolean networks [32].
Although multistate systems can be converted into Boolean networks [10], the conversion adds new artificial nodes, and thus models might lose their intuitiveness or become more computationally expensive to analyze. In this manuscript, we have presented a method for extending mixed-state networks into a multistate system in a natural way that preserves steady states. Namely, we convert the system to a PDS over a finite field. We also present control methods based on computational algebra that can generate new steady states, destroy existing steady states, or avoid regions in the state space. Importantly, our control strategies allows for the targeting of both nodes and edges of the wiring diagram.
We also note that, in theory, computing the Gröbner basis for a system of polynomial equations can be computationally expensive (with doubly exponential complexity). However, for many biological systems, computing the Gröbner basis can be achieved in a reasonable time [18] and the computational cost does not seem to correlate with the size of the network but with the average connectivity [47].
Although the methods presented here focus on steady states and synchronous updating schedules many biological systems are modeled with asynchronous schedules and stochastic methods. Developing efficient algebraic methods for the control of asynchronous and stochastic multiscale models presents a rich opportunity for the development of methods more widely applicable to biological systems. We will explore computational algebra methods for general updating schedules in the future.

Conclusion
In many settings in biology and, in particular, in biomedicine, the ultimate goal of an investigation is the solution of a control problem, and mathematical modeling can be a helpful tool in this endeavor. Modeling dynamic biological networks using systems of ordinary differential equations has the advantage that the modeler has ready access to well-developed rigorous mathematical control theory tools. These work well, when enough quantitative information is available. Some control problems, however, such as the ones considered here, are of a more qualitative nature, such as modifying a network as to change its steady state structure in specified ways. These are not as well suited to a control approach based on differential equations, but fit naturally into a discrete modeling framework. There are now several rigorous approaches to control for discrete models that collectively allow a rigorous solution of a range of control problems. This paper adds another methodology to this collection, based on the principle of representing networks through a collection of polynomials over a finite field, which makes available the algorithms, software, and mathematical tools of computational algebra and algebraic geometry for the solution of a wide range of related problems. The algorithm in this paper makes this approach available for general (deterministic) discrete dynamic networks.

Supplementary Materials
All scripts used to generate the data in this manuscript can be found in the first author's github repository https://github.com/luissv7/multistatepdscontrol. The software used for computations of Gröbner bases was a Macaulay2 Docker Image v: 1.14 [15].

Appendix
In this appendix we denote finite fields with either F q or F p , where p is assumed to be a prime number while q is assumed to be a power of a prime number.

Converting mixed-state models into polynomial dynamical systems
Let q be the smallest number which is a power of a prime number such that q ≥ |X i | for all i. Consider the finite field F = F q .
We can identify X i → F by an injective map ι i for i from 1 to n. Let ι = (ι 1 , . . . , ι n ). We can now consider the dynamical system F as a subsystem of a dynamical system F : F n → F n as follows.
Notice that α i essentially "crushes" the points in F − ι i (X i ) into a constant in ι i (X i ).
Here ι = (ι 1 , ι 2 ), α = (α 1 , α 2 ), and F 2 Veliz-Cuba et al. [48] previously used a similar transformation for a finite field of prime order, F p , where the elements outside of F p − ι(X) were sent into the "largest" element (p − 1). However, in a general finite field F q , there is no adequate concept of the "largest element". Notice thatF(x 1 , . . . , x n ) = (x 1 , . . . , x n ) if and only if (x 1 , . . . , x n ) is in the image of ι and (ι −1 1 (x 1 ), . . . , ι −1 n (x n )) is a fixed point of F. In particular, we can now "extend" the discrete dynamical system F to a discrete dynamical system F : F n → F n without changing the dynamics of the original system. 7.2 An approach for deriving a polynomial dynamical system from a mixed-state dynamical system A common approach to representing mixed-state dynamical systems is to give Boolean expressions for when a certain node will attain a given value based on the state of the other nodes [56,37]. For example, in the signaling network model presented in [37], the rule for representing how E2F3 attains values 1 or 2 are shown in Table 4.
In the case that some of the variables are Boolean (can only take one of two values), and the other variables are in a set of the same prime cardinality q, we can convert to a polynomial dynamical system over F q . If a variable x i was Boolean to start with, we replace x i with x q−1 i . For a variable, x i that was not Boolean, we can write the polynomial representation by taking advantage of indicators functions q j (x) = (Π i∈F q ,i = j (x − i)) q−1 for j ∈ F q . For example, if a variable appears in a Boolean expression as x i = j, then we substitute that variable with (Π i∈F q ,i = j (x − i)) q−1 . Recall that the operator AND is equivalent to the product over F 2 , the operator OR is equivalent to the operator (x, y) → x + y − (x + y) and NOT is equivalent to x → 1 + x. Over F q , we define x AND y to be (x, y) → (x · y) q−1 , NOT x to be x → 1 − x q−1 and x OR y to be (x, y) → −(x · y) q−1 + x q−1 + y q−1 .

Continuity condition and steady states
The continuity condition is a restriction that the state of each variable does not change by more than one unit at each time step (see e.g. [3] for details). Intuitively, the continuity condition represents that a biological quantity cannot suddenly go from high to low (or low to high) without reaching an intermediate step.
Here we show that the continuity condition on polynomial dynamical systems used in [3] does not change steady states. Fix a prime p and consider the finite field k = F p . Fix the notation x = (x 1 , . . . , x i−1 , x i , x i+1 , . . . , x n ), and let F i := f i (x). We will always assume that the representative for x i is in the set {0, 1, · · · , p − 1}.
We will say that f : k[x 1 , · · · , x n ] → k n is continuous if Any PDS F : k n → k n can be made continuous by consideringF : k n → k n wherê F i = h • (F i × π i ) where π i is the projection onto the ith coordinate.
Theorem 1 Let F : k n → k n be a polynomial dynamical system over a finite field k and letF : k n → k n be the polynomial dynamical system where the continuity condition has been applied to F. Then the set of fixed points of F, FIX(F) is equal to FIX(F) Proof Let x ∈FIX(F), π i : k n → k be the projection onto the i − th coordinate.
Now, if x ∈ FIX(F), we have h(F i (x), x i ) = x i for all i. This can only happen if As a result, we have FIX(F) =FIX(F).