Quantifying the Total Effect of Edge Interventions in Discrete Multistate Networks

TL;DR

This paper develops polynomial formulas to quantify the global impact of edge interventions in discrete multistate networks, especially Boolean networks, using canalizing functions, and validates these formulas with random network simulations.

Contribution

It introduces a polynomial normal form for nested canalizing functions and formulas to count transition changes due to edge deletions, advancing the analysis of intervention effects in network models.

Findings

Formulas accurately estimate maximum transition changes upon edge deletion.

Canalizing structure influences the number of affected transitions.

Upper bounds are validated and characterized in random networks.

Abstract

Developing efficient computational methods to assess the impact of external interventions on the dynamics of a network model is an important problem in systems biology. This paper focuses on quantifying the global changes that result from the application of an intervention to produce a desired effect, which we define as the total effect of the intervention. The type of mathematical models that we will consider are discrete dynamical systems which include the widely used Boolean networks and their generalizations. The potential interventions can be represented by a set of nodes and edges that can be manipulated to produce a desired effect on the system. We use a class of regulatory rules called nested canalizing functions that frequently appear in published models and were inspired by the concept of canalization in evolutionary biology. In this paper, we provide a polynomial normal form based on the canalizing properties of regulatory functions. Using this polynomial normal form, we give a set of formulas for counting the maximum number of transitions that will change in the state space upon an edge deletion in the wiring diagram. These formulas rely on the canalizing structure of the target function since the number of changed transitions depends on the canalizing layer that includes the input to be deleted. We also present computations on random networks to compare the exact number of changes with the upper bounds provided by our formulas. Finally, we provide statistics on the sharpness of these upper bounds in random networks.

Equations

$$\mathbf{F}=(f_1,\dots ,f_n):X\to X$$

$${\mathcal{F}}_j(x,u)=f_j(x_{j_1},\dots ,(1-Q_{i,j}(u))x_i,\dots ,x_{j_m}).$$

$${\mathcal{F}}_j(x,u)=f_j(x_{j_1},\dots ,x_i=0,\dots ,x_{j_m}).$$

$${\mathcal{F}}_j(x,u)=f_j(x_{j_1},\dots ,x_i,\dots ,x_{j_m}).$$

$$\begin{array}[]{l}\mathcal{F}_{j}(x,u)=\\ f_{j}(x_{j_{1}},\dots,(1-Q_{i,j}(u))x_{i}+aQ_{i,j}(u),\dots,x_{j_{m}}).\end{array}$$

$${\mathcal{F}}_j(x,u)=f_j(x_{j_1},\dots ,x_i=a,\dots ,x_{j_m}).$$

$${\mathcal{F}}_j(x,u)=f_j(x_{j_1},\dots ,x_i,\dots ,x_{j_m}).$$

$$\begin{array}[]{l}f(x_{1},\dots,x_{n})=\\ \\ \left\{\begin{array}[]{l}b_{1},\text{ if}\ x_{\sigma(1)}\in S_{1},\\ b_{2},\text{ if}\ x_{\sigma(1)}\notin S_{1},x_{\sigma(2)}\in S_{2},\\ \vdots\\ b_{k},\text{ if}\ x_{\sigma(1)}\notin S_{1},\dots,x_{\sigma(k)}\in S_{k},\\ g\neq b_{k},\text{ if}\ x_{\sigma(1)}\notin S_{1},\dots,x_{\sigma(k)}\notin S_{k},\end{array}\right.\end{array}$$

$$f(x_1,x_2)=1+2x_1^2+2x_2+2x_1^2{x}_2+2x_2^2.$$

$$\widetilde{Q}_{S}(x)=\biggl{\{}\begin{array}[]{ll}1&\text{if}\ x\notin S,\\ 0&\text{if}\ x\in S.\end{array}$$

$$\begin{array}[]{l}f(x_{1},\dots,x_{n})=M_{1}(M_{2}(\dots(M_{r-1}(M_{r}P_{C}+\\ \\ \quad\quad\quad B_{r})+B_{r-1})\dots)+B_{2})+B_{1},\\ \end{array}$$

$$g=M_2(\dots (M_{r-1}(M_r{P}_C+B_r)+B_{r-1})\dots )+B_2.$$

$$\begin{array}[]{l}f(x_{1},\dots,x_{n})=M_{1}(M_{2}(\dots(M_{r-1}(B_{r+1}M_{r}+\\ \quad B_{r})+B_{r-1})\dots)+B_{2})+B_{1},\end{array}$$

$$\begin{array}[]{l}f(x_{1},x_{2},x_{3},x_{4})=1+x_{1}^{2}+x_{1}^{2}x_{2}+2x_{1}^{2}x_{2}^{2}+\\ x_{1}^{2}x_{2}x_{3}+2x_{1}^{2}x_{2}^{2}x_{3}+x_{1}^{2}x_{2}x_{4}+2x_{1}^{2}x_{2}^{2}x_{4}.\end{array}$$

$$f=M_1(M_2(P_C+1)+1)+1),$$

$$\begin{array}[]{l}f(x_{1},x_{2},x_{3},x_{4})=\\ \\ \left\{\begin{array}[]{l}1,\text{ if}\ x_{1}\in S_{1}=\{0\},\\ 2,\text{ if}\ x_{1}\notin S_{1},x_{2}\in S_{2}=\{0,1\},\\ P_{C},\text{ if}\ x_{1}\notin S_{1},x_{2}\notin S_{2}.\end{array}\right.\end{array}$$

$$\begin{array}[]{l}f_{t}(x_{1},\dots,x_{n})=M_{1}^{t}(M_{2}^{t}(\dots(M_{r-1}^{t}(M_{r}^{t}P_{C}+\\ \\ \quad\quad\quad\quad B_{r})+B_{r-1})\dots)+B_{2})+B_{1},\end{array}$$

$$\begin{array}[]{l}p^{n-k_{1}-\cdots-k_{\ell}}\left(p-\left|S\right|\right)^{k_{1}+\cdots+k_{\ell}}.\end{array}$$

$$\begin{array}[]{l}p^{n-k_{1}-\cdots-k_{\ell}}\cdot\left(\displaystyle\prod_{i=1}^{\ell-1}\leavevmode\nobreak\ \leavevmode\nobreak\ \prod_{j=1}^{k_{i}}(p-|S_{j,t}|)\right)\\ \left(\displaystyle\prod_{\begin{subarray}{c}j=1\\ j\neq s\end{subarray}}^{k_{\ell}}(p-|S_{j,t}|)\right)(p-R),\end{array}$$

$$R=\left\{\begin{array}[]{cc}|S_{s,t}|&\textrm{if $0\in S_{s,t}$}\\ p-|S_{s,t}|&\textrm{if $0\notin S_{s,t}$.}\end{array}\right.$$

$$\begin{array}[]{l}\frac{p^{n}}{p^{k_{1}+\cdots+k_{\ell-1}}}\cdot\left(\displaystyle\prod_{i=1}^{\ell-1}\leavevmode\nobreak\ \leavevmode\nobreak\ \prod_{j=1}^{k_{i}}(p-|S_{j,t}|)\right)\\ \left(\frac{1}{p^{k_{\ell}-1}}\displaystyle\prod_{\begin{subarray}{c}j=1\\ j\neq s\end{subarray}}^{k_{\ell}}(p-|S_{j,t}|)\right)\frac{1}{p}(p-R)=\\ p^{n-k_{1}-\cdots-k_{\ell}}\cdot\left(\displaystyle\prod_{i=1}^{\ell-1}\leavevmode\nobreak\ \leavevmode\nobreak\ \prod_{j=1}^{k_{i}}(p-|S_{j,t}|)\right)\\ \left(\displaystyle\prod_{\begin{subarray}{c}j=1\\ j\neq s\end{subarray}}^{k_{\ell}}(p-|S_{j,t}|)\right)(p-R),\end{array}$$

$$R=\left\{\begin{array}[]{cc}|S_{s,t}|&\textrm{if $a\in S_{s,t}$}\\ p-|S_{s,t}|&\textrm{if $a\notin S_{s,t}$.}\end{array}\right.$$

$$p^{n-k_1-\cdots -k_r-1}\left(\prod _{i=1}^r\text{leavevmode}\prod _{j=1}^{k_i}(p-|S_{j,t}|)\right)(p-1).$$

$$p^{n-k_1-\cdots -k_r-1}\left(\prod _{i=1}^r\text{leavevmode}\prod _{j=1}^{k_i}(p-|S_{j,t}|)\right)(p-R)$$

$$R=\left\{\begin{array}[]{cc}|S_{s,t}|&\textrm{if $0\in S_{s,t}$}\\ 1&\textrm{if $0\notin S_{s,t}$.}\end{array}\right.$$

$$\begin{array}[]{l}p^{n-k_{1}-\cdots-k_{r}-2}\left(\displaystyle\prod_{i=1}^{r}\leavevmode\nobreak\ \prod_{j=1}^{k_{i}}(p-|S_{j,t}|)\right)\\ (p-|S_{s,t}|)(p-1).\end{array}$$

$$\begin{array}[]{l}p^{n-k_{1}-\cdots-k_{r}-c-1}\left(\displaystyle\prod_{i=1}^{r}\leavevmode\nobreak\ \prod_{j=1}^{k_{i}}(p-|S_{j,t}|)\right)\\ \prod_{i=1}^{c}(p-|S_{s_{i},t}|)(p-1),\end{array}$$

Full text

aDepartment of Mathematics, University of Kentucky,

Lexington, KY 40506-0027, USA

bDepartment of Mathematics, California Polytechnic

State University, San Luis Obispo, CA 93407-0403, USA

1 Introduction

Boolean networks (BN) have been proposed as an appropriate framework for modeling the state of cells due to their simplicity and the variety of tools available for model analysis Veliz-Cuba et al. (2010); Zañudo & Albert (2015). However, some complex gene interactions cannot be represented in the Boolean setting and several generalizations of the Boolean approach have been developed Thomas & D’Ari (1990). Multistate models, a generalization of the BN framework, where the genes can attain more than two states have been proposed as appropriate models for capturing complex gene expression patterns, such as consideration of three states (low, medium, and high). We note that while in theory it is possible to develop models where the variables can take on any number of possible states (possibly only restricted by the requirement that it needs to be a power of a prime number so that the domain can have the structure of a finite field), in practical applications this number is typically small – rarely larger than 5.

A Gene Regulatory Network (GRN) is a representation of the intricate relationships among genes, proteins, and other compounds that are responsible for the expression levels of mRNA and proteins. Boolean networks have been successfully used to model and study the properties of GRN Albert & Othmer (2003); Li et al. (2004). In particular, Boolean canalizing rules were introduced by S. Kauffman and collaborators Kauffman et al. (2003, 2004) and reflect the concept of canalization in evolutionary biology that Waddington pioneered in 1942 Waddington (1957) – that organisms evolve developmental robustness, producing an invariant phenotype even under genetic or environmental perturbations.

In this article, we study the network-wide effect of an experimental intervention that either prevents a regulation from happening or silences a node. Such intervention is modeled through edge deletion and can be achieved via therapeutic drugs that target a specific gene interaction Choi et al. (2012); Campbell & Albert (2019). In Murrugarra & Dimitrova (2015) we introduced methods for quantifying side effects in Boolean networks. However, many of the more recently published discrete dynamical models include variables that take on more than two states due to the need for capturing mechanisms that are not binary in nature Zañudo, Scaltriti & Albert (2017); Remy et al. (2015); Chifman et al. (2017); Espinosa-Soto et al. (2004); Thieffry & Thomas (1995). Consequently, Boolean nested and partially nested canalizing functions were generalized to multistate Murrugarra & Laubenbacher (2011, 2012); Kadelka et al. (2017) which enables the possibility of capturing more complex interactions among the genes in the network. Such functions can be viewed as a discrete dynamical system with a stratified structure which consists of hierarchical layers of variables according to their relative influence over the system dynamics.

The ability to quantify the global changes in the dynamics of the network after an external perturbation has important applications. In the presence of external network modifications where the topology of the network changes but the attractor structure remains unchanged, it is still desirable to quantify the changes in other aspects of the dynamics such as the transient time. For instance, in evolutionary biology to simulate evolution one often evolves an ensemble of networks (by performing mutations, crossover, selection, etc.) for many generations Wagner (1996). At the end of the simulations, one measures the changes in the evolved networks to compare with the features of the original ensemble. In Wagner (1996), after simulated evolution, the evolved networks had similar features to the original ones such as the number of attractors and basin sizes. One feature that had changed is the transient time Wagner (1996). The theoretical tools presented in this paper will be useful to measure global changes even if the attractor structure is preserved after an intervention.

There are several published control methods for Boolean networks such as Stable Motifs Zañudo & Albert (2015), Feedback Vertex Sets Zañudo, Yang & Albert (2017), Minimal Hitting Sets Vera-Licona et al. (2013); Klamt et al. (2006), and several others Qiu et al. (2014); Li et al. (2015); Poret & Boissel (2014); Gates & Rocha (2016); Zañudo, Scaltriti & Albert (2017); Sordo Vieira et al. (2019). While these control methods focus on finding control interventions that satisfy a control objective (e.g., to drive the system into a specific attractor), there are very few studies of the total extent of the consequences of applying a certain control action (beyond satisfying the control objective). This paper contributes methods for measuring the impact of the control actions on the dynamics of multistate networks. The type of theoretical tools presented here can help to discriminate control actions with minimal effect on the state space. That is, even if we have different control candidates that can achieve a certain objective, they might have different impact on the dynamics of the network and we might be interested in distinguishing the control action that produces the least changes in the dynamics of the network.

The rest of the paper is structured as follows. In Section 2, we introduce discrete dynamical systems and their representation as polynomial dynamical systems. In Section 3 we define the control actions for multistate networks. In Section 4 we provide a polynomial normal form for discrete functions and then we use this representation to derive a set of formulas for counting the maximum number of transitions in the state space upon edge deletions. In Section 5 we present computational results for random networks. Finally, in Section 6, we provide the conclusions of the paper.

2 Background

A discrete dynamical system can be defined as a dynamical system that is discrete in time as well as in variable states. More formally, consider a collection $x_{1},\ldots,x_{n}$ of variables, each of which can take on values in finite sets $X_{1},\dots,X_{n}$. Let $X=X_{1}\times\dots\times X_{n}$ be their Cartesian product. A discrete dynamical system in the variables $x_{1},\ldots,x_{n}$ is a function

[TABLE]

where each coordinate function $f_{i}$ is a discrete function on a subset of $\{x_{1},\dots,x_{n}\}$ which represents how the future value of the $i$-th variable depends on the present values of the variables. If $X_{i}=\{0,1\}$, then each $f_{i}$ is a Boolean rule and $\mathbf{F}$ is a Boolean network.

In this article, for the purpose of exploiting the algebraic properties of discrete functions, it is assumed that the variables $x_{1},\ldots,x_{n}$ take on values from a finite field $\mathbb{F}$. Then using the fact that any discrete function $f_{i}:\mathbb{F}^{n}\rightarrow\mathbb{F}$ can be represented as a polynomial in $x_{1},\dots,x_{n}$, that is $f_{i}\in\mathbb{F}[x_{1},\dots,x_{n}]$, the discrete network can be represented as $\mathbf{F}=(f_{1},\dots,f_{n}):\mathbb{F}^{n}\rightarrow\mathbb{F}^{n}$ where each $f_{i}\in\mathbb{F}[x_{1},\dots,x_{n}]$. If any of the variables $x_{1},\ldots,x_{n}$ take on values from a set that cannot be directly identified with a finite field, then it is straightforward to embed the system $\mathbf{F}:X\rightarrow X$ into a system $\hat{\mathbf{F}}:\mathbb{F}^{n}\rightarrow\mathbb{F}^{n}$, where $X\subset\mathbb{F}^{n}$, while preserving the attractor structure of $\mathbf{F}$; see Veliz-Cuba et al. (2010).

Given a discrete network $\mathbf{F}=(f_{1},\dots,f_{n})$, a directed graph $\mathcal{W}$ on $n$ nodes $x_{1},\ldots,x_{n}$ is associated to $\mathbf{F}$ as follows: there is a directed edge in $\mathcal{W}$ from $x_{j}$ to $x_{i}$ if $x_{j}$ appears in $f_{i}$, i.e. $x_{j}$ is in the support of $f_{i}$, written $x_{j}\in supp(f_{i})$. In the context of a molecular network model, this graph represents the wiring diagram of the network.

The dynamics of a discrete network is given by the difference equation $x(t+1)=\mathbf{F}(x(t))$; that is, the dynamics is generated by iteration of $\mathbf{F}$. More precisely, the dynamics of $\mathbf{F}$ is represented by the state space graph $S$, defined as the graph with vertices in $\mathbb{F}^{n}$ which has an edge from $x\in\mathbb{F}^{n}$ to $y\in\mathbb{F}^{n}$ if and only if $y=\mathbf{F}(x)$. In this context, the problem of finding the states $x\in\mathbb{F}^{n}$ where the system will get stabilized is of particular importance. The collection of these special points of the state space are called attractors of a discrete network and elements of the attractors may include steady states (fixed points), where $\mathbf{F}(x)=x$, or cycles, where $\mathbf{F}^{r}(x)=x$ for some integer $r>1$. Attractors in network modeling might represent cell types Kauffman (1969) or cellular states such as apoptosis, proliferation, or cell senescence Huang (1999); Shmulevich & Dougherty (2010).

3 Methods

Network interventions can be modeled through edge and node manipulations and can be achieved via therapeutic drugs that target a specific gene interaction Choi et al. (2012); Campbell & Albert (2019). In Murrugarra & Dimitrova (2015); Murrugarra et al. (2016) we provided definitions for these actions in Boolean networks. These definitions are usually used for encoding the control parameters with the purpose of identifying control targets as shown in Murrugarra et al. (2016). In this paper we will consider the deletion and constant expression of edges in the multistate setting.

Throughout the paper, $S_{i,j}$ will be a subset of $\mathbb{F}$ and $Q_{i,j}(u)$ will be the indicator functions of $S_{i,j}$. That is, they return 1 when $u$ is in the set and 0 when $u$ is not. The index $i$ indicates the node $x_{i}$ from which the edge begins and the second index $j$ is used when necessary to mark the function under consideration.

3.1 Edge Control in Multistate Networks

In the Boolean setting, the deletion of an edge 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.

Definition 3.1** (Edge Deletion).**

Consider the edge $x_{i}\rightarrow x_{j}$ in a wiring diagram. For $u\in S_{i,j}$, the control of the edge $x_{i}\rightarrow x_{j}$ consists of manipulating the input variable $x_{i}$ for $f_{j}$ in the following way:

[TABLE]

For each value of $u\in\mathbb{F}$ we have the following control settings:

• •

For $u\in S_{i,j}$,

[TABLE]

That is, the control is active and the action represents the removal of the edge $x_{i}\rightarrow x_{j}$.

• •

For $u\notin S_{i,j}$,

[TABLE]

That is, the control is not active.

We will also consider the constant expression of edges, which we define as follows.

Definition 3.2** (Constant expression).**

Consider the edge $x_{i}\rightarrow x_{j}$ in a wiring diagram and $a\in\mathbb{F}$. For $u\in S_{i,j}$, the control of the edge $x_{i}\rightarrow x_{j}$ consists of manipulating the input variable $x_{i}$ for $f_{j}$ in the following way:

[TABLE]

For each value of $u\in\mathbb{F}$ we have the following settings:

• •

For $u\in S_{i,j}$,

[TABLE]

That is, the control is active and the action represents the constant expression (to $a$) of the edge $x_{i}\rightarrow x_{j}$.

• •

For $u\notin S_{i,j}$,

[TABLE]

That is, the control is not active.

Remark 3.3**.**

Node control of $x_{i}$ can be defined as setting the function $f_{i}$ to a constant $a\in\mathbb{F}$.

4 Results

In this section we present a definition of $k$-canalizing functions for the multistate case and then we characterize this functions in terms of layers of canalizations. Subsequently, we use this canalizing layers representation to derive an upper bound for the number of changes in the state space of a discrete system upon an edge deletion in the wiring diagram.

4.1 Multistate $k$-Canalizing Functions

In the following definition, we assume that $\sigma$ is a permutation on $\{1,\dots,n\}$.

Definition 4.1**.**

The function $f:\mathbb{F}^{n}\rightarrow\mathbb{F}$ is a $k$-canalizing function in the variable order $x_{\sigma(1)},\dots,x_{\sigma(k)}$ with canalizing input sets $S_{1},\dots,S_{k}\subset\mathbb{F}$ and canalizing output values $b_{1},\dots,b_{k}\in\mathbb{F}$ if it can be represented in the form

[TABLE]

where $g=g(x_{\sigma(k+1)},\dots,x_{\sigma(n)})$ is a multistate function on $n-k$ variables. When $g$ is not a canalizing function, the integer $k$ is the canalizing depth of $f$. If $g$ is not a constant function, then $g$ is called the core function of $f$ and is denoted by $P_{C}$.

Remark 4.2**.**

Note that in Definition 4.1 we require that the function $g$ be unique when all the canalizing variables are not in their corresponding canalizing input sets. As a result, a function could be canalizing but not $1$-canalizing, see Example 4.3.

Example 4.3**.**

Let $\mathbb{F}=\{0,1,2\}$ and $n=2$. Consider the function

[TABLE]

For this function $x_{2}$ is canalizing (with $S_{1}=\{2\}$) because $f(x_{1},2)=1$. However, $f$ is not a $1$-canalizing function because $f(x_{1},0)=1+2x_{1}^{2}\neq 2+x_{1}^{2}=f(x_{1},1)$. Thus, even though $x_{2}$ is canalizing for $f$ , the function $f$ has no layers of canalization. Thus, $P_{C}=f$.

4.2 Layers of canalization in multistate networks

In Theorem 4.4 we provide a polynomial normal description of discrete functions. Basically, this theorem gives a partition of the inputs of the function into canalizing and non-canalizing variables and, within the canalizing ones, we categorize the input variables into layers of canalization. This theorem is a generalization of a theorem in He & Macauley (2016) from Boolean to the multistate case.

Let $S\subset\mathbb{F}$ be a subset of $\mathbb{F}$ and $\widetilde{Q}_{S}(u)$ be the indicator function of the complement of $S$. That is,

[TABLE]

Theorem 4.4**.**

Every multistate function can be uniquely written as

[TABLE]

where $M_{i}=\displaystyle\prod_{j=1}^{k_{i}}\widetilde{Q}_{S_{i,j}}$, $d=k_{1}+\cdots+k_{r}$ is the canalizing depth, $P_{C}$ is a polynomial that has no canalizing variables, $B_{1},B_{2},\dots,B_{r}\in\mathbb{F}$, and $B_{r}\neq 0$. Each variable $x_{i}$ appears in exactly one of the $M_{1},M_{2},\dots,M_{r},P_{C}$.

Proof.

If $f(x_{1},\dots,x_{n})$ is non-canalizing, then $P_{C}=f$. If $f(x_{1},\dots,x_{n})$ is canalizing, then we proceed by induction. For $n=1$, if $f$ is canalizing in $x_{i}$ but not $1$-canalizing in $x_{i}$, then we set $P_{C}=f$. If $f$ is $1$-canalizing in $x_{i}$, then it can be written as $f=\widetilde{Q}_{S_{1}}(x_{i})+B_{1}$ for some set $S_{1}\subset\mathbb{F}$. Then $f$ has the form of Equation 2 by setting $M_{1}=\widetilde{Q}_{S_{1}}(x_{i})$ and $P_{C}=1$. For $n=2$, if $f(x_{i},x_{j})$ is not $1$-canalizing on any of its variables, then we set $P_{C}=f$. If $f$ is $1$-canalizing on $x_{i}$, then $f$ can be written as $f(x_{i},x_{j})=M_{1}(x_{i})g(x_{j})+B_{1}$ for some $g(x_{j})$. Then $f$ has the form of Equation 2 by setting $P_{C}=g$. Now assume that Equation 2 is true for any canalizing function that is essential in at most $n-1$ variables (that is, for all functions that depend in at most $n-1$ variables). Let $f$ be a function that is essential in $n$ variables. If $f$ is not $1$-canalizing on any of its variables, then we set $P_{C}=f$. If $f$ is $1$-canalizing in $x_{i_{1}},\dots,x_{i_{k_{1}}}$, then $f=M_{1}g+B_{1}$, where $M_{1}$ is the product of indicator functions of the complements of sets $S_{i_{1}},\dots,S_{i_{k_{1}}}\subset\mathbb{F}$ and $g$ has $n-k_{1}$ variables. If $g$ has no canalizing variables, then $f$ has the form of Equation 2 with $P_{C}=g$. If $g$ is canalizing, then by the inductive hypothesis $g$ can be written as

[TABLE]

Thus, $f$ has the form of Equation 2. ∎

Remark 4.5**.**

For a multistate nested canalizing function, the formula in Equation 2 reduces to

[TABLE]

as was shown in Kadelka et al. (2017).

In the following example we describe a $2$-canalizing function with noncanalizing variables.

Example 4.6**.**

Let $\mathbb{F}=\{0,1,2\}$ and $n=4$. Consider the function

[TABLE]

The function $f$ can be written as in Equation 2 as

[TABLE]

where $M_{1}=\widetilde{Q}_{S_{1}}(x_{1})=x_{1}^{2}$, $S_{1}=\{0\}$, $M_{2}=\widetilde{Q}_{S_{2}}(x_{2})=x_{2}+2x_{2}^{2}$, $S_{2}=\{0,1\}$, and $P_{C}=x_{3}+x_{4}$. Thus $f$ has two layers and two noncanalizing variables. Note that $f$ can also be written as in Equation 1 as

[TABLE]

4.3 Upper bounds

Using the polynomial normal form of multistate functions in Theorem 4.4, we derive a set of formulas for counting the maximum number of transitions that will change in the state space upon an edge deletion in the wiring diagram. The formulas presented here are generalizations from the Boolean case to the multistate setting of the formulas we presented in Murrugarra & Dimitrova (2015).

For the next theorem, we are going to assume that the functions of the discrete network $\mathbf{F}=(f_{1},\ldots,f_{n}):\mathbb{F}^{n}\to\mathbb{F}^{n}$ are written in the format of Theorem 4.4. That is, for $t=1,\dots,n$ the coordinate function $f_{t}$ has the following form,

[TABLE]

where $M^{t}_{i}=\displaystyle\prod_{j=1}^{k_{i}}\widetilde{Q}_{S_{j,t}}$, $d=k_{1}+\cdots+k_{r}$ is the canalizing depth, $P_{C}$ is a polynomial with no canalizing variables, $B_{1},B_{2},\dots,B_{r}\in\mathbb{F}$, and $B_{r}\neq 0$. Each variable $x_{i}$ appears in exactly one of $M_{1}^{t},M_{2}^{t},\dots,M_{r}^{t},P_{C}$.

Remark 4.7**.**

Note that the function $f_{t}$ has $r$ layers and there are $k_{i}$ variables in each layer for $i=1,\dots,r$.

In the following theorem, we assume that the canalizing input sets are all of the same size for all the variables. In Theorem 4.10 we study the general case where the canalizing input sets of the variables can be different.

Theorem 4.8**.**

Let $\mathbf{F}=(f_{1},\ldots,f_{n}):\mathbb{F}^{n}\to\mathbb{F}^{n}$ be a multistate network where $f_{t}$ is a $k$-canalizing function written as in Eq. 4 with $k_{1},\ldots,k_{r}$ the numbers of variables in layers $1,\ldots,r$, respectively. Let $x_{s}$ be in the $\ell^{th}$ layer, where $\ell\leq r$ and $r$ is the number of layers. Suppose that all canalizing variables have the same canalizing input set $S$ and that $0\in S$. Then, the maximum number of transitions in the state space that will change upon deletion of $x_{s}\to x_{t}$ is given by

[TABLE]

Proof.

Let $m=k_{1}+\cdots+k_{\ell}$. The number of input vectors where the other canalizing variables (not $x_{s}$) of $f_{t}$ do not take on their canalizing input is $\left(p-\left|S\right|\right)^{m-1}$. For these input vectors, if $x_{s}$ was already set to 0 or to any other of its canalizing values in $S$, then the output of $f_{t}$ will not change as a result of deleting $x_{s}\to x_{t}$. Finally, since we have $n-m$ non-canalizing variables, the total number of input vectors for which the output of $f_{t}$ can possibly change is $\left(p-\left|S\right|\right)^{m-1}(p-\left|S\right|)p^{n-m}=\left(p-\left|S\right|\right)^{m}p^{n-m}$. ∎

Remark 4.9**.**

Note that from Equation 5 that the number of variables in each layer affects the number of changes and that there are potentially more changes when the deletion happens in a more dominant layer, see examples 5.1-5.2.

Theorem 4.10**.**

Let $\mathbf{F}=(f_{1},\ldots,f_{n}):\mathbb{F}^{n}\to\mathbb{F}^{n}$ be a multistate network where $f_{t}$ is a $k$-canalizing function written as in Eq. 4 with $k_{1},\ldots,k_{r}$ the numbers of variables in layers $1,\ldots,r$, respectively. Let $x_{s}$ be in the $\ell^{th}$ layer, $\ell\leq r$ and $r$ is the number of layers. The maximum number of transitions in the state space that will change upon deletion of $x_{s}\to x_{t}$ is given by

[TABLE]

where

[TABLE]

Proof.

The strategy is to first count the number of inputs that do not contain values from the canalizing sets of the variables in the first $\ell-1$ layers (that do not contain $x_{s}$). Thus, the term in the first line of Equation 7 counts the number of non-canalizing inputs in the previous layers to the layer containing $x_{s}$; the term inside the second set of parentheses of Equation 7 counts the number of non-canalizing inputs of the variables (except of $x_{s}$) in the layer containing $x_{s}$; the last term in Equation 7 counts the number of non-canalizing inputs of $x_{s}$. For the last term, notice that deleting $x_{s}\to x_{t}$ results in setting $x_{s}=0$ in $f_{t}$. If 0 is in the canalizing set of $x_{s}$, $S_{s,t}$, then the rest of the values in $S_{s,t}$ will yield the same output as 0. Since $|S_{s,t}|/p$ of the input values in the transition table of $f_{t}$ contain a canalizing value for $x_{s}$, it is the remaining $\frac{p-|S_{s,t}|}{p}$ of the table that can potentially change as a result of the edge deletion. On the other hand, if $0\notin S_{s,t}$, then it is the inputs not in $S_{s,t}$ that have the potential to change the output as a result of deleting $x_{s}\to x_{t}$ which constitutes $1/p$ of the transition table, with $\frac{p-1}{p}$ of the table that can potentially change as a result of the edge deletion. Thus, to obtain Equation 6 we multiply the following expressions:

[TABLE]

∎

Remark 4.11**.**

The bound in Equation 6 is sharp. 2. 2.

When $p=2$, the formula in Equation 6 reduces to $2^{n-k_{1}-k_{2}-\cdots-k_{r}}$. 3. 3.

If instead of edge deletion, we consider constant expression to $a\in\mathbb{F}$ (see Section LABEL:subsec:ce) of $x_{s}\to x_{t}$, then the formula in Equation 6 remains the same except for $R$ which becomes

[TABLE]

Proposition 4.12**.**

Let $\mathbf{F}=(f_{1},\ldots,f_{n}):\mathbb{F}^{n}\to\mathbb{F}^{n}$ be a multistate network where $f_{t}$ is written as in Equation 4. Let $x_{s}\in supp(P_{C})$. The maximum number of transitions in the state space that will change upon deletion of $x_{s}\to x_{t}$ is

[TABLE]

Remark 4.13**.**

This upper bound is sharp. 2. 2.

When $p=2$, the expression reduces to $2^{n-d-1}$, where $d$ is the canalizing depth. 3. 3.

If $f$ has no canalizing variables, then the formula in Equation 8 reduces to $p^{n-1}(p-1)$.

Proposition 4.14**.**

Let $\mathbf{F}=(f_{1},\ldots,f_{n}):\mathbb{F}^{n}\to\mathbb{F}^{n}$ be a multistate network where $f_{t}$ is written as in Equation 4. If $P_{C}$ is canalizing but not $1$-canalizing with canalizing variable $x_{s}$ and input set $S_{s,t}$, then there are two cases to consider:

The deletion of $x_{s}\to x_{t}$ will result in up to

[TABLE]

transitions, where

[TABLE] 2. 2.

Let $x_{a}\in supp(P_{C})$ that is not canalizing. Then the maximum number of transitions in the state space that will change upon deletion of $x_{a}\to x_{t}$ is

[TABLE]

Remark 4.15**.**

This upper bound is sharp. 2. 2.

If $P_{C}$ has more than one canalizing variables (but it is still not 1-canalizing), then the formula in Equation 10 becomes

[TABLE]

where each $x_{s_{i}}$ is a canalizing variable and $c$ is the number of canalizing variables of $P_{C}$.

5 Applications

To provide further insights into the results presented above and to illustrate the use of the formulas here we present numerical results for random networks where we compare the exact number of changes to the upper bounds provided by the formulas.

For the next examples, we generated random networks with scale-free structure using the Barabasi-Albert algorithm Barabasi & Albert (1999). We note that the Barabasi-Albert algorithm produces undirected edges but for our examples we need directed edges. Thus, for each undirected edge between $x_{i}$ and $x_{j}$, we converted the edge $x_{i}-x_{j}$ into either $x_{i}\rightarrow x_{j}$ or $x_{j}\rightarrow x_{i}$ at random.

Example 5.1** (Boolean Case).**

In order to calculate the exact number of changes we use random networks with 10 nodes. For each network, we selected the node with the maximum in-degree and generated a random partition of its inputs to assign the canalizing layers. The functions of the other nodes were generated at random.

In Figure 1 we show statistics for the number of changes in the first layer. The average maximum in-degree for the networks in Figure 1 is $4.14$ ($std=1.07$). The average number of variables in the first layer is $2.61$ ($std=1.5$). The average number of changes in the first layer is $221.28$ ($std=168.363$) and the average upper bound is $259.04$ ($std=201.562$).

In Figure 2 we present statistics of the number of changes as well the difference between the exact number of changes and the upper bound provided by the formulas. For these networks, in about 40% of the cases the upper bounds match the exact number of changes.

In Figure 3 we show statistics for the number of changes in the second layer. The average maximum in-degree for the networks in Figure 3 is $4.8$ ($std=1.07$). The average number of variables in the first layer is $1.67$ ($std=0.93$). The average number of variables in the second layer is $2.13$ ($std=0.75$). The average number of changes in the second layer is $86.0$ ($std=62.7$) and the average upper bound is $100.64$ ($std=78.88$).

In Figure 4 we present statistics of the number of changes as well the difference between the exact number of changes and the upper bound provided by the formulas. For these networks, in about 50% of the cases the upper bounds match the exact number of changes.

From Figures 2 and 4, we see that the bounds are slightly more accurate when the edge intervention happens in a less dominant layer.

Example 5.2** (Multistate Case).**

Here we also use random networks with scale-free structure with $p=3$ and $n=10$ nodes. For each network, we selected the node with the maximum in-degree and generated a random partition of its inputs to assign the canalizing layers. The functions of the other nodes were generated at random.

In Figure 5 we show statistics for the number of changes in the first layer. The average maximum in-degree for the networks in Figure 5 is $4.05$ ($std=0.88$). The average number of variables in the first layer is $2.28$ ($std=1.16$). The average number of changes in the first layer is $17303.2$ ($std=14739.4$) and the average upper bound is $17792.5$ ($std=15220.6$).

In Figure 6 we present statistics of the difference between the upper bounds provided by the formulas and the exact number of changes. For these networks, in about 75% of the cases the upper bounds match the exact number of changes.

In Figure 7 we show statistics for the number of changes in the second layer. The average maximum in-degree for the networks in Figure 7 is $4.84$ ($std=0.94$). The average number of variables in the first layer is $1.7$ ($std=0.86$). The average number of variables in the second layer is $1.93$ ($std=0.97$). The average number of changes in the second layer is $3946.68$ ($std=3789.71$) and the average upper bound is $4056.48$ ($std=3969.7$).

In Figure 8 we present statistics of the difference between the the upper bound provided by the formulas and the exact number of changes. For these networks, in about 80% of the cases the upper bounds match the exact number of changes.

From Figures 6 and 8, we see that the bounds are slightly more accurate when the edge intervention happens in a less dominant layer.

6 Conclusions

In this paper we present practical methods for quantifying the global changes that result from an application of an external intervention in the network, which we called the total effect of the intervention. We emphasized that, while there are several methods for identifying control targets in discrete networks, there have been very few studies focusing on the total extent of the consequences of applying a certain control action (beyond satisfying the control objective). This paper contributes methods for measuring the number of changed transitions in the state space upon the application of an edge control in multistate networks. The approach is based on a polynomial normal form description of discrete functions that provides a way of categorizing the inputs of the function and therefore of quantifying their impact on the dynamics of the network. The main computational cost of our approach is in obtaining the canalizing layers format of the functions which we used to derive our formulas. Once the functions are written in the layers format, it is straightforward to apply the formulas to calculate the upper bound. We note that obtaining the layers format could be a formidable task with exponential complexity on the number of inputs in the worst case but for the type of networks we are studying (i.e. biological networks) we believe that our approach is still practical. We applied our methods to randomly generated multistate models and verified that in many cases the upper bounds provided by our formulas were accurate. We also observed that the upper bounds tend to be more accurate when the edge interventions happen in a less dominant layer.