Solving Equations on Discrete Dynamical Systems (Extended version)

TL;DR

This paper introduces a practical method for solving certain decidable polynomial equations over discrete dynamical systems, enabling hypothesis verification in biological models with promising scalability.

Contribution

It proposes a novel approach to solve specific decidable equations in discrete dynamical systems, facilitating hypothesis testing in biological modeling.

Findings

Approach effectively reduces complex equations to simpler forms.

Experimental results show good scalability of the method.

Method enables enumeration of all solutions for hypothesis verification.

Abstract

Boolean automata networks, genetic regulation networks, and metabolic networks are just a few examples of biological modelling by discrete dynamical systems (DDS). A major issue in modelling is the verification of the model against the experimental data or inducing the model under uncertainties in the data. Equipping finite discrete dynamical systems with an algebraic structure of commutative semiring provides a suitable context for hypothesis verification on the dynamics of DDS. Indeed, hypothesis on the systems can be translated into polynomial equations over DDS. Solutions to these equations provide the validation to the initial hypothesis. Unfortunately, finding solutions to general equations over DDS is undecidable. In this article, we want to push the envelope further by proposing a practical approach for some decidable cases in a suitable configuration that we call the Hypothesis Checking. We demonstrate that for many decidable equations all boils down to a "simpler" equation. However, the problem is not to decide if the simple equation has a solution, but to enumerate all the solutions in order to verify the hypothesis on the real and undecidable systems. We evaluate experimentally our approach and show that it has good scalability properties.

Figures

Tables

$x_1$	$x_2$	…	$x_p$	$x_1$	$x_2$	…	$x_p$	…	$x_1$	$x_2$	…	$x_p$
$y_1$	$y_2$	…	$y_{q_1}$	$y_q$	$y_1$	…	$y_{q-2}$	…	…	…	…	$y_q$

$x_1$	$x_2$	…	$x_i$	$x_{i+1}$	…	$x_p$
$y_1$	$y_2$	…	$y_q$	$y_1$	…	$y_q$

Table 3. Table 1 : Final data-structure storing all the decompositions, each solution for each value and at each step, the set of all solutions for a given value.

Node	Splits	Node solution	Subtree solutions set
6	[3,3][2,2,2]	$C_6^1$	$\{C_6^1,C_3^2,C_1^1\oplus C_2^1\oplus C_3^1,C_3^1\oplus C_1^3,$ $C_2^1\oplus C_1^4,C_1^6,C_2^3,C_1^2\oplus C_2^2\}$
3	[2,1]	$C_3^1$	$\{C_3^1,C_1^1\oplus C_2^1,C_1^3\}$
2	[1,1]	$C_2^1$	$\{C_1^2,C_2^1\}$
1	$\varnothing$	$C_1^1$	$\left\{C_1^1\right\}$

Table 4. Table 2 : Final data-structure storing all the decomposition, each solution for each value and at each step, the set of all solutions for a given value.

Node	Splits	Node solution	Subtree solutions set
5	[4,1]	$\left\{\right\}$	$\left\{\right\}$
4	[2,2]	$\left\{\right\}$	$\left\{C_4^2\right\}$
2	[1,1]	$C_4^1$	$\left\{C_4^1\right\}$
1	$\varnothing$	$\left\{\right\}$	$\left\{\right\}$

Table 5. Table 3 : Final data-structure storing all the decomposition, each solution for each value and at each step, the set of all solutions for a given value.

Node	Splits	Node solution	Subtree solutions set
12	[6,6]	$\left\{\right\}$	$\{C_3^4\oplus C_6^4,C_3^{12},C_6^6,C_3^6\oplus C_6^3,$ $C_6^2\oplus C_3^8,C_3^2\oplus C_6^5,C_6^1\oplus C_3^{10}\}$
6	[3,3] [2,2,2]	$\left\{\right\}$	$\{C_3^6,C_6^2\oplus C_3^2,C_3^4\oplus C_6^1,C_6^3\}$
3	[2,1]	$\left\{\right\}$	$\{C_3^3,C_6^1\oplus C_3^1\}$
2	[1,1]	$C_6^1$	$\{C_6^1,C_3^2\}$
1	$\varnothing$	$C_3^1$	$\left\{C_3^1\right\}$

Equations

$$a_1\cdot X_1+a_2\cdot X_2+\dots +a_k\cdot X_k=C$$

$$P(x_1,\dots ,x_n)=Q(x_1,\dots ,x_n).$$

$$C_p^m\odot C_q^n=C_{\mathrm{lcm}(p,q)}^{m\times n\times \gcd (p,q)}.$$

$$C_p^1\odot C_q^1=C_{\mathrm{lcm}(p,q)}^{\gcd (p,q)}$$

$$C_p^{mn}=n\cdot C_p^m.$$

$${\left(S\right)}^n=\underset{i=1}{\overset{m}{⨁}}{C}_{p_i}^{p_i^{n-1}}\oplus \underset{\begin{array}{c}{k}_1+k_2+...+k_m=n\\ 0\le k_1,k_2,...,k_m<n\end{array}}{⨁}\left(\genfrac{}{}{0ex}{}{n}{k_1,k_2,...,k_m}\right)C_l^{g\cdot {\prod }_{\begin{array}{c}t=1\\ k_t\ne 0\end{array}}^m{p}_t^{k_t-1}}.$$

$$(S{)}^n=(C_{p_1}^1\oplus C_{p_2}^1\oplus ...\oplus C_{p_m}^1{)}^n=\underset{k_1+k_2+...+k_m=n}{⨁}\left(\genfrac{}{}{0ex}{}{n}{k_1,k_2,...,k_m}\right)\underset{t=1}{\overset{m}{⨀}}(C_{p_t}^1{)}^{k_t}=$$

$$=\underset{i=1}{\overset{m}{⨁}}(C_{p_i}^1{)}^n\oplus \underset{\begin{array}{c}{k}_1+k_2+...+k_m=n\\ 0\le k_1,k_2,...,k_m<n\end{array}}{⨁}(\genfrac{}{}{0ex}{}{n}{k_1,k_2,...,k_m})\underset{t=1}{\overset{m}{⨀}}(C_{p_t}^1{)}^{k_t}$$

$$\underset{i=1}{\overset{m}{⨁}}(C_{p_i}^1{)}^n\oplus \underset{\begin{array}{c}{k}_1+k_2+...+k_m=n\\ 0\le k_1,k_2,...,k_m<n\end{array}}{⨁}(\genfrac{}{}{0ex}{}{n}{k_1,k_2,...,k_m})\underset{t=1}{\overset{m}{⨀}}(C_{p_t}^1{)}^{k_t}=$$

$$=\underset{i=1}{\overset{m}{⨁}}(C_{p_i}^1{)}^n\oplus \underset{\begin{array}{c}{k}_1+k_2+...+k_m=n\\ 0\le k_1,k_2,...,k_m<n\end{array}}{⨁}(\genfrac{}{}{0ex}{}{n}{k_1,k_2,...,k_m})C_{l(p_1,p_2,...,p_m,k_1,k_2,...,k_m)}^{g(p_1,p_2,...,p_m,k_1,k_2,...,k_m)\cdot {\prod }_{t=1}^m{p}_t^{k_t-1}}=$$

$$=\underset{i=1}{\overset{m}{⨁}}{C}_{p_i}^{p_i^{n-1}}\oplus \underset{\begin{array}{c}{k}_1+k_2+...+k_m=n\\ 0\le k_1,k_2,...,k_m<n\end{array}}{⨁}\left(\genfrac{}{}{0ex}{}{n}{k_1,k_2,...,k_m}\right)C_{l(p_1,p_2,...,p_m,k_1,k_2,...,k_m)}^{g(p_1,p_2,...,p_m,k_1,k_2,...,k_m)\cdot {\prod }_{t=1}^m{p}_t^{k_t-1}}.$$

$$(\underset{j=1}{\overset{S_1}{⨁}}{C}_{p_{1j}}^1\odot X_1)\oplus (\underset{j=1}{\overset{S_2}{⨁}}{C}_{p_{2j}}^1\odot X_2)\oplus \dots \oplus (\underset{j=1}{\overset{S_k}{⨁}}{C}_{p_{kj}}^1\odot X_k)=\underset{j=1}{\overset{m}{⨁}}{C}_{q_j}^{n_j}$$

$$C_p^1\odot X=C_q^n.$$

Full text

Solving Equations on Discrete Dynamical Systems

Alberto Dennunzio

  

Enrico Formenti

  

Luciano Margara

  

Valentin Montmirail

  

Sara Riva

1 Scientific Background

Boolean automata networks have been heavily used in the study of systems biology [2, 6]. The main drawback of the approach by automata network is in the very first step, namely when one induces the network from the experiments. Indeed, most of the time the knowledge of the network is partial and hypotheses are made about its real structure. Those hypotheses must be verified either by further experiments or by the study of the dynamical evolution of the network compared to the expected behaviour provided by the experimental evidences.

In [3], an abstract algebraic setting for representing the dynamical evolution of finite discrete dynamical systems has been proposed. In the following, we denote by $R$, the commutative semi-ring of the DDS.

The basic idea is to identify a discrete dynamical system with the graph of its dynamics (finite graphs having out-degree exactly $1$) and then define operations $+$ and $\cdot$ which compose dynamical systems to obtain larger ones.

Indeed, a discrete dynamical system (DDS) is a structure $\left\langle\chi,f\right\rangle$ where $\chi$ is a finite set called the set of states and $f:\chi\to\chi$ is a function called the next state map. Any DDS $\left\langle\chi,f\right\rangle$ can be identified with its dynamics graph which is a structure $G\equiv\left\langle V,E\right\rangle$ where $V=\chi$ and $E=\left\{(a,b)\in V\times V,\,f(a)=b\right\}$. From now on, when speaking of a DDS, we will always refer to its dynamics graph.

Given two DDS $G_{1}=\left\langle V_{1},E_{1}\right\rangle$ and $G_{2}=\left\langle V_{2},E_{2}\right\rangle$ their sum $G_{1}+G_{2}$ is defined as $\left\langle V_{1}\cupdot V_{2},E_{1}\cupdot E_{2}\right\rangle$, where $\cupdot$ denotes the disjoint union. The product $G_{1}\cdot G_{2}$ is the structure $\left\langle V^{\prime},E^{\prime}\right\rangle$ where $V^{\prime}=V_{1}\times V_{2}$ and $E^{\prime}=\{((a,x),(b,y))\in V^{\prime}\times V^{\prime},(a,b)\in E_{1}\text{ and }(x,y)\in E_{2}\}$. It is easy to see that $F\equiv\left\langle\chi,+,\cdot\right\rangle$ is a commutative semiring in which $\left\langle\emptyset,\emptyset\right\rangle$ is the neutral element w.r.t. $+$ and $\left\langle\left\{a\right\},\left\{(a,a)\right\}\right\rangle$ is the neutral element w.r.t. $\cdot$ operation.

Now, consider the semiring $R[X_{1},X_{2},\ldots,X_{n}]$ of polynomials over $R$ in the variables $X_{i}$, naturally induced by $R$. Let us go back to our initial motivation. Assume that some parts of the overall dynamics $a_{1},a_{2},\ldots,a_{k}$ are known, then the following equation represents a hypothesis on the overall structure of the expected dynamical system $C$ on the basis of the known data $a_{1},\ldots,a_{k}$, where all the coefficients, variables and $C$ are DDS.

[TABLE]

The hypotheses are verified whenever the previous equation admits a solution, therefore providing a way to solve such equation can be used to check hypotheses against a given discrete dynamical system. For the sake of clarity, we denote our unknown variables as $X_{i}$, whereas they, in fact, represent any monomial of the form $x_{i}^{w_{i}}$. The following fundamental result states that solving polynomial equations over DDS is not an easy task.

Theorem 1.1** ([3]).**

Given two polynomials $P(x_{1},\ldots,x_{n})$ and $Q(x_{1},\ldots,x_{n})$ over $R[x_{1},\ldots,x_{n}]$, consider the following equation

[TABLE]

The problem of finding a solution to Equation 2 is undecidable. Moreover, if Equation 2 is linear or quadratic, then finding a solution is in \NP. Finally, when $P(x)=const$, where the polynomial is in a single variable and all its coefficients are systems consisting of self-loops only, the equation is solvable in polynomial time.

According to Theorem 1.1, solving polynomial equations of the type $P(x_{1},\ldots,x_{n})=const$ is in \NP even for quadratic polynomials. In order to overcome this issue, one can follow at least two strategies: either further constrain the polynomials or solve approximated equations which can provide information on the real solution.

In this article, we follow the second option. Indeed, we focus on strongly connected components (SCC) of the dynamics graph. Recall that SCC represents a very important feature in finite DDS since they are the attracting sets. These sets contain the asymptotic information about system evolution.

2 Methods

In the dynamics graph, each component of a system can be divided in two parts: the transient part and the periodic part, see [5] for more details. A point $x\in X$ of a discrete dynamical system $\langle X,f\rangle$ belongs to a cycle if there exists a positive number $p\in\mathbb{N}$ such that $f^{p}(x)=x$. The smallest $p$ is the period of the cycle, and $x$ is periodic. The periodic part is the set of nodes periodic. All the others nodes are transient, but in this work, X is a finite set hence any state $x$ is ultimately periodic and in each component of the graph there is only one cycle of length at least 1.

Every finite DDS can be described as a sum of single components, and every component can be described, for our purposes, with the length of its period (strongly connected components in dynamics graphs are cycles). The transient part of a component is not relevant for the result of the sum and product operations when the equation is over SCC.

A single component of period $p$ is denoted $C^{1}_{p}$, while $C^{n}_{p}$ means that there are $n$ components of period p in the system. Therefore, if a system is composed by $n$ components, each of period $p_{i}$ with $i\in\left\{1,\ldots,n\right\}$, then $\bigoplus\limits_{i=1}^{n}C^{1}_{p_{i}}$ completely describes the system where $\bigoplus$ denotes the sum of components since each component has only one period (see Figure 1).

Remark 2.1**.**

When a system has several components with the same period, then their representation can be added. As an example, we have $C^{1}_{2}\oplus C^{1}_{2}=C^{2}_{2}$. Otherwise, the sum $\oplus$ consists of a concatenation of components.

3 Contributions

From now on, $\bar{R}$ will indicate the restriction of $R$ to systems made by strongly connected components only. First, we need to adapt the definition of product between two DDS in terms of components and their period.

Proposition 3.1**.**

Consider a system composed by $m$ components of period $p$, namely $C^{m}_{p}$, multiplied by a system with $n$ components of period $q$, namely $C^{n}_{q}$, the result of the product operation depends only on the length of the periods of the components involved according to the following formula, with $m,n\in\mathbb{N}$ and $m,n\geq 1$

[TABLE]

Proof 3.2**.**

Given two discrete dynamical systems $\langle X,f\rangle$ and $\langle Y,g\rangle$, where the first system has only one component of period $p$ and the second has only one component of period $q$, let us prove that:

[TABLE]

We know that a product operation corresponds to a Cartesian product between $X$ and $Y$ (Given two discrete dynamical systems $\langle X,f\rangle$ and $\langle Y,g\rangle$, their product is the dynamical system $\langle X\times Y,f\times g\rangle$ where $\forall(x,y)\in X\times Y,(f\times g)(x,y)=(f(x),g(y)))$. There are two possible cases:

• •

$gcd(p,q)=1$.

In this case the larger period of the two is not able to represent the smallest cyclic behavior inside it, consequently it obtains a single period containing all the Cartesian product.

• •

$gcd(p,q)\neq 1$, with $p>q$.

In this case the cycles of period $p$ and $q$ arrive at one point to be represented by a cycle of length $lcm(p,q)$ but this means that the elements of this cycle are only a subset of the Cartesian product. For this reason $gcd(p,q)={{p\cdot q}\over{lcm(p,q)}}$ components are generated.

In the case of $m$ or $n$ different from $1$, this means that each product operation is done for each of these components, so in general the result is duplicated $m\cdot n$ times.

One can also simplify the parameter of a component. The following definition provides a formula to compact the notation of a DDS with $n$ identical components.

Definition 3.3**.**

Consider a single component $C_{p}^{m}$, then $\forall n,m,p\in\mathbb{N}\setminus\left\{0\right\}$ it holds

[TABLE]

Let us remind that each $X_{i}$ represents, in fact, a variable $x_{i}^{w_{i}}$. Therefore, it is necessary to know how we can retrieve the solutions for the original $x_{i}$. To do so, we will use the following lemma:

Proposition 3.4**.**

Given a system composed by $m$ components of period $p_{i}\in\mathbb{N}$, with $p_{i}>0$ for all $i\in\{1,...,m\}$, let $g(p_{1},p_{2},...,p_{m},k_{1},k_{2},...,k_{m})$ be the gcd between the $p_{i}$ for which $k_{i}\not=0$ and let $l(p_{1},p_{2},...,p_{m},k_{1},k_{2},...,k_{m})$ be the lcm between the $p_{i}$ for which $k_{i}\not=0$. Consider a system $S\equiv C_{p_{1}}^{1}\oplus C_{p_{2}}^{1}\oplus...\oplus C^{1}_{p_{m}}$. Then,

[TABLE]

Proof 3.5**.**

Using the multinomial theorem one finds

[TABLE]

[TABLE]

The resulting Formula 5 is obtained by extrapolating the cases in which a $k_{i}=n$. Another transformation is possible according to Proposition 3.1.

[TABLE]

[TABLE]

[TABLE]

For $k$ equal to [math] we assume that $(S)^{0}$ is equal to $C^{1}_{1}$, the neutral element of the product operation. Let us go back to Equation 1 which is the problem that we want to solve. It can be rewritten as follows:

[TABLE]

with $S_{i}$, the number of different components in the system $i$, $p_{ij}$ is the value of the period of the $j^{th}$ component in the system $i$. In the right term, there are $m$ different periods, where for the $j^{th}$ different period, $n_{j}$ is the number of components, and $q_{j}$ the value of the period. However, Equation (6) is still hard to solve. We can simplify it performing a contraction step which consists in cutting Equation (6) into two simpler equations: $(C_{p_{11}}^{1}\odot X_{1})=W$, where $W\subseteq\bigoplus\limits_{i=1}^{m}C_{q_{i}}^{n_{i}}$ and $((C_{1}^{1}\odot Y)=\bigoplus\limits_{i=1}^{m}C_{q_{i}}^{n_{i}}\setminus W)$ with $Y=(\bigoplus\limits_{i=2}^{S_{1}}C_{p_{1i}}^{1}\odot X_{1})\oplus(\bigoplus\limits_{j=1}^{S_{2}}C_{p_{2j}}^{1}\odot X_{2})\oplus\ldots\oplus(\bigoplus\limits_{j=1}^{S_{k}}C_{p_{kj}}^{1}\odot X_{k})$. By applying recursively a contraction step on all the partitions of $W$ and on the second equation obtained (i.e. the one containing $Y$) one finds that, solving Equation (6) boils down to solving multiple times the following type of equation:

[TABLE]

If the variable $X$ presents a power different from one, it is possible use the Lemma 3.4 in order to study the squared by the power.

However, equations of the shape of Equation 7 will be numerous therefore an efficient practical algorithm able to enumerate all its solutions is needed. In fact, we can propose the following bounds to know how many times equations of the shape Equation 7 are solved with the following lemma:

Proposition 3.6**.**

Let us denote by $Z$ the number of times that we will solve equation of the shape Equation 7, we have the following: $\prod_{i=1}^{m}\binom{n_{i}+\sum_{j=1}^{k}S_{j}-1}{\sum_{j=1}^{k}S_{j}-1}\cdot m\leq Z\leq\prod_{i=1}^{m}\binom{n_{i}+\sum_{j=1}^{k}S_{j}-1}{\sum_{j=1}^{k}S_{j}-1}\cdot m\cdot\sum_{j=1}^{k}S_{j}$.

The intuition is as follows: the contraction step is necessary to study all the possible ways to produce the right term with the components in the left part of the equation. Accordingly, it is necessary to understand the number of possible decompositions of the right term to discover the bounds for the number of the executions of the colored-tree method (a decomposition corresponds to assign a subset of the components of the right part to a product operation between a variable and a known component). For each period a Star and Bars decomposition is applied (we redirect the reader unfamiliar with the Star and Bars decomposition to [4]).

Proof 3.7**.**

In general for a fixed $q_{i}$, the $n_{i}$ components are divided into $\sum_{j=1}^{k}S_{j}$ groups, in this case there are $\binom{n_{i}+\sum_{j=1}^{k}S_{j}-1}{\sum_{j=1}^{k}S_{j}-1}$ different ways for dividing the components. Therefore, we can rewrite the lemma as follows: $m\leq\frac{Z}{\prod_{i=1}^{m}\binom{n_{i}+\sum_{j=1}^{k}S_{j}-1}{\sum_{j=1}^{k}S_{j}-1}}\leq m\cdot\sum_{j=1}^{k}S_{j}$. And now, toward a contradiction for the lower bound. Let us assume that we can solve less than $m$ equations. This implies that we solve less equations than the number of different periods on the right term. Contradiction, we need at least all of them (not necessary all their combinations) to determine the solution of the equation. And now, toward a contradiction again to prove the upper-bound. Firstly, we know for all the components in the right term there are $\prod_{i=1}^{m}\binom{n_{i}+\sum_{j=1}^{k}S_{j}-1}{\sum_{j=1}^{k}S_{j}-1}$ feasible divisions. Now, let us assume that in the worst case, for each coefficient the product operation must produce more than one components of each possible periods in the right term. This is a contradiction from the definition of the equation, where all the components must all have a different period. The second possibility to go beyond this bounds is that it would exists more $S_{i}$ than the one present in the equation, again a contradiction by definition of the equation. Therefore, we know that we have: $\prod_{i=1}^{m}\binom{n_{i}+\sum_{j=1}^{k}S_{j}-1}{\sum_{j=1}^{k}S_{j}-1}\cdot m\leq Z\leq\prod_{i=1}^{m}\binom{n_{i}+\sum_{j=1}^{k}S_{j}-1}{\sum_{j=1}^{k}S_{j}-1}\cdot m\cdot\sum_{j=1}^{k}S_{j}$, for $Z$ being the number of times that we will solved equation of the shape Equation 7.

4 The Colored-Table Method

First of all, let us formally define the problem and analyze its complexity.

Definition 4.1** (\DSECP).**

The (finite) Discrete Dynamical Systems Solving Equations on Components Problem is a problem which takes in input $C^{1}_{p}$ and $C^{n}_{q}$ and outputs the list of all the solutions $X$ to the equation $C^{1}_{p}\odot X=C^{n}_{q}$.

Solving \DSECP is hard but still tractable. Indeed, the following lemma classifies our problem in \EnumP. Recall that \EnumP is the complexity class of enumeration problems for which a solution can be verified in polynomial time [7]. It can be seen as the enumeration counterpart of the \NP complexity class.

Lemma 4.2**.**

\DSECP

is in \EnumP.

Proof 4.3**.**

One just needs to be able to check if a given value is a solution in polynomial time. This can be done in linear time using Lemma 3.1.

Notation.

For any $n,p,q\in\mathbb{N}^{\star}$, let $T_{p,q}^{n}$ denote the set of solutions of Equation (7) and $S_{p,q}^{n}$ the set of solutions returned by the colored-tree method.

The colored-tree method is pretty involved, we prefer start to illustrate it by an example.

Example 4.4**.**

Consider the following equation $C_{6}^{1}\odot X=C_{6}^{6}$. The algorithm consists in two distinct phases: tree building and solution aggregation. In the first phase, the algorithm enumerates all the divisors $\mathcal{D}$ of $6$ i.e. $\left\{6,3,2,1\right\}$. It then applies a making-change decomposition algorithm (MCDA) [1] in which the total sum is $6$ and the allowed set of coins is $\mathcal{D^{\prime}}=\mathcal{D}\setminus\left\{6\right\}$. MCDA decomposes $6$ as $3+3$ (which is an optimal decomposition). MCDA is then applied recursively (always using $\mathcal{D}\setminus\left\{i\right\}$ as the set of coins to decompose $i$). We obtain $(6=3+3)$, $(3=2+1)$ and $(2=1+1)$ as reported in Table 1.

At this point, a check is performed to ensure that all possible ways of decomposing $6$ using $\mathcal{D^{\prime}}$ are present in the tree. In our case, we already have $[3,3]$ found by the first run of MCDA. We also found: $[3,2,1]$, $[2,2,1,1]$, $[1,1,2,1,1]$, $[1,1,1,1,1,1]$ by the recursive application of MCDA. By performing the check, we discover that the decomposition of $6$ as $[2,2,2]$ is not represented in the current tree. For this reason, $[2,2,2]$ is added to the set of decompositions of $6$ as illustrated in Figure 2, it is assigned a new color and a recursive application of MCDA is started on the newly added nodes. A new check ensures that all decompositions are present. This ends the building phase. The resulting tree is reported in Figure 2.

*After this first phase of construction of the tree, the aggregation of solutions starts. Remark that each node $m$ represents the equation $C^{1}_{p}\odot X=C^{m}_{q}$ that we call the node equation. The single component solution is called the node solution and it is obtained thanks to Lemma 3.1, $C^{1}_{{\frac{q}{p}}\times m}$ whenever a feasible solution exists i.e. if $\gcd(p,{\frac{q}{p}}\times m)=m$ and $\operatorname{lcm}(p,{\frac{q}{p}}\times m)=q$. For example, for $m=3$ one finds $x=C_{3}^{1}$. To find all the solutions for the current node one must also take the Cartesian product of the solutions sets in the subtrees of the same color and then the union of the solution sets of nodes of different colors (different splits). All the solutions can be found in Table 1. *

Example 4.5**.**

Consider the equation $C_{2}^{1}\odot X=C_{4}^{5}$. In the first phase, the algorithm enumerates all the divisors $\mathcal{D}$ of $4$ i.e. $\left\{4,2,1\right\}$. It then applies a making-change decomposition algorithm (MCDA) [1]. MCDA decomposes $5$ as $4+1$ (which is an optimal decomposition). MCDA is then applied recursively always using $\mathcal{D}\setminus\left\{i\right\}$ as the set of coins to decompose $i$. We obtain $(5=4+1)$, $(4=2+2)$ and $(2=1+1)$ as reported in Table 2.

At this point, a check is performed to ensure that all possible ways of decomposing $5$ using $\mathcal{D}\setminus\left\{i\right\}$ as the set of coins to decompose $i$. In our case, we already have $[4,1]$ found by the first run of MCDA. We also found: $[2,2,1]$, $[2,1,1,1]$, $[1,1,1,1,1]$ by the recursive application of MCDA. By performing the check, we discover that all the possible decompositions of $5$ are represented in the current tree. This ends the building phase. The resulting tree is reported in Figure 3.

After this first phase of construction of the tree, the aggregation of solutions starts. In this case the tree presents only one color. Remark that if in the cartesian product a empty set is involved, the result of the operation is the empty set. For example, for $m=2$ , one has that the node solution is $C_{4}^{1}$. From the subtrees of the node one finds a empty set, but with the union of the solution of the node, the subtree solutions set for $m=2$ is $\left\{C_{4}^{1}\right\}$. Moreover, the final solution set for the node $5$ is the empty set, in fact in the Cartesian product $m=1$ is involved (empty set). In this case the method return a empty set of solutions, that represents the impossibility of the equation.

Example 4.6**.**

Consider the equation $C_{2}^{1}\odot X=C_{6}^{12}$. In the first phase, the algorithm enumerates all the divisors $\mathcal{D}$ of $6$ i.e. $\left\{6,3,2,1\right\}$. It then applies a making-change decomposition algorithm (MCDA) [1]. MCDA decomposes $12$ as $6+6$ (which is an optimal decomposition). MCDA is then applied recursively always using $\mathcal{D}\setminus\left\{i\right\}$ as the set of coins to decompose $i$. We obtain $(12=6+6)$, $(6=3+3)$, $(3=2+1)$ and $(2=1+1)$ as reported in Table 3.

At this point, a check is performed to ensure that all possible ways of decomposing $12$ using $\mathcal{D^{\prime}}$ is present in the tree. In our case, the decomposition of $6$ in $[2,2,2]$ is added in "each occurrence" of $6$. This ends the building phase. The resulting tree is reported in Figure 4.

After this first phase of construction of the tree, the aggregation of solutions starts. To find the solutions for the current node one must also take the Cartesian product of the solutions sets in the subtrees of the same color and then the union of the solution sets of nodes of different colors (different splits). For example, for $m=12$ (i.e. the root node), the cartesian product between $6$ and $6$ is computed, but for $m=6$ (in each occurrence) two cartesian operations and a union are necessary. Therefore, the final solution set for the node $12$ is $\left\{C^{4}_{3}\oplus C^{4}_{6},C^{12}_{3},C^{6}_{6},C^{6}_{3}\oplus C^{3}_{6},C^{2}_{6}\oplus C^{8}_{3},C^{2}_{3}\oplus C^{5}_{6},C^{1}_{6}\oplus C^{10}_{3}\right\}$.

Although we can describe our algorithm with a pseudocode, and then we can sketch some proofs about its soundness, completeness and termination.

The Lisiting LABEL:coloredtable presents the procedure using some particular functions:

• •

generateNewNodes adds the elements of the split, the node necessary in order to decompose but not yet represented as nodes in the nodes set.

• •

MCDA computes the optimal solutions of the making-change problem for a node value and a set of coins.

• •

computeSingleSolution returns the node solution for a node equation represented with a node.

• •

checkRepresented check if all the possible decomositions of the root are represented, otherwise add the corrisponding sub-tree.

• •

IncreaseOrder permutes the row of the table in the increasing order according to the value of the nodes.

Now we can sketch some proofs about its soundness, completeness and termination.

Proposition 4.7** (Soundness).**

For all $n,p,q\in\mathbb{N}^{\star},\;S_{p,q}^{n}\subseteq T_{p,q}^{n}$.

Proof 4.8**.**

Let us prove the soundness by induction on the depth of the tree from leaves to root. Induction base: if there is only one step, we know by Lemma 3.1, that a solution found is feasible iff $gcd(p,{\frac{q}{p}}\times m)=m$ and $lcm(p,{\frac{q}{p}}\times m)=q$, and because there is only one leaf in the base, we therefore, obtain all the solutions. Induction hypothesis: let us assume that we have all the possible solutions at a depth $n$ and let us show that we can obtain all the solutions at a depth $n+1$. Induction step: It is easy to see that a solution exists if and only if it comes from a decomposition. Thus, by performing a Cartesian product between the set of solutions at depth $n$ (which is true by IH) and the node solution (which is true by Induction base, since the node can be seen as a leaf), we know that we will obtain all the solution coming from the possible decomposition in the sub-tree. If a solution is coming from another sub-tree, since we perform an exhaustive check where we assign a different color to the other sub-tree, we know again, by IH and because we are taking the union of all the possible solutions, that we have all the possible solutions at a depth $n+1$.

Proposition 4.9** (Completeness).**

For all $n,p,q\in\mathbb{N}^{\star},\;T_{p,q}^{n}\subseteq S_{p,q}^{n}$.

Proof 4.10**.**

By contradiction, let us assume that there exists a solution $r\in T_{p,q}^{n}$ and that $r\not\in S_{p,q}^{n}$. This means that the colored-tree method does not return it. This implies that it exists a decomposition of $n$, which leads to $r$, such that this decomposition is not in the tree. This is impossible since, an exhaustive check is performed to assure that all the decompositions are there. Therefore, all solutions are returned.

Proposition 4.11** (Termination).**

The colored-tree method always terminates.

Proof 4.12**.**

The building phase always terminates since the colored-tree has maximal depth $\mathcal{D^{\prime}}=div(q,n)$ and the number of different possible colors is bounded by $2^{k}$ where $k$ is the size of the multi-set containing $n/p_{i}$ copies of the divisor $p_{i}$ per each divisor in $\mathcal{D^{\prime}}$. The aggregation phase always terminates since it performs a finite number of operations per each node of the colored tree.

Now that we have defined the problem, its complexity and a sound and complete algorithm to solve it. It is time to experimentally evaluate it in order to study its scalability.

5 Experimental Evaluations

The colored-tree method provides a complete set of solutions of simple equations of type Equation 7. Its complexity can be experimentally measured counting the number of nodes in the colored tree.

Figure 5 shows how the complexity grows as a function of $n$ and $q$. For this case, we set $p=q$ to ensure that we always have at least one solution and therefore a tree-decomposition. Notice that, in some cases, the complexity is particularly high due to specific analytical relations between the input parameters that we are going to study in the future. Notice also that our method seems to have a weakness when $q$ is an even number. This is easily explained: in many cases, all the divisors can be expressed by the other ones. Therefore the check that ensures that all the decompositions are present is particularly time- and memory-consuming.

Since there is no other competitor algorithm at the best of our knowledge, we compared the colored-tree method to a brute force algorithm. We test our algorithm with $n$ from $1$ to $20$, $p$ is also from $1$ to $20$ and at any time, $p=q$. Results are reported in Figure 6. As expected, the colored-tree method outperforms the brute force solution, sometimes with many orders of magnitude faster. However, when the input equation has small coefficients, the colored-tree method performs worse. This can be explained considering that building the needed data structures requires a longer time than the execution of the brute force algorithm.

6 Conclusion

Questions about boolean automata networks, used in biological modelling for genetic regulatory networks and metabolic networks, can be rewritten as equations over DDS using the formalism introduced in Dennunzio et al. in [3]. They argued that polynomial equations are a convenient tool for the analysis of the dynamics of a system. However, algorithmically solving such equations is an unfeasible task. In this article, we propose a practical way to partially overcome those difficulties using a couple of approaches which aims at studying separately the number of component (i.e. the number of attractors) and the length of their periods. This paper proposes an algorithm for the number of components of the solution of a polynomial equation over finite DDS.

One of the core routines of the algorithm uses a brute force check for the make-change problem which clearly affects the overall performances. Therefore, a natural research direction consists in finding a better performing routine. One possibility would consider parallelisation since a large part of the computations are strictly indipendent. Another interesting research direction consists inbetter understanding the computational complexity of the \DSECP. We are still working to improve the performances of the algorithm to have stronger scalability properties in the perspective of providing a handy tool which can be exploited by bioinformaticians to actually solve the Hypothesis Checking problem in their context.

Acknowledgments

This work has been partly funded by IDEX UCA${}^{\textsc{jedi}}$.