 Research
 Open Access
 Published:
Reducing Boolean networks with backward equivalence
BMC Bioinformatics volume 24, Article number: 212 (2023)
Abstract
Background
Boolean Networks (BNs) are a popular dynamical model in biology where the state of each component is represented by a variable taking binary values that express, for instance, activation/deactivation or high/low concentrations. Unfortunately, these models suffer from the state space explosion, i.e., there are exponentially many states in the number of BN variables, which hampers their analysis.
Results
We present Boolean Backward Equivalence (BBE), a novel reduction technique for BNs which collapses system variables that, if initialized with same value, maintain matching values in all states. A largescale validation on 86 models from two online model repositories reveals that BBE is effective, since it is able to reduce more than 90% of the models. Furthermore, on such models we also show that BBE brings notable analysis speedups, both in terms of state space generation and steadystate analysis. In several cases, BBE allowed the analysis of models that were originally intractable due to the complexity. On two selected case studies, we show how one can tune the reduction power of BBE using modelspecific information to preserve all dynamics of interest, and selectively exclude behavior that does not have biological relevance.
Conclusions
BBE complements existing reduction methods, preserving properties that other reduction methods fail to reproduce, and vice versa. BBE drops all and only the dynamics, including attractors, originating from states where BBEequivalent variables have been initialized with different activation values The remaining part of the dynamics is preserved exactly, including the length of the preserved attractors, and their reachability from given initial conditions, without adding any spurious behaviours. Given that BBE is a modeltomodel reduction technique, it can be combined with further reduction methods for BNs.
Background
Boolean networks (BNs) are a popular model in systems biology where the dynamics is qualitatively associated with two levels. These may express, for instance, on/off behavior in gene regulation or high/low concentrations of molecular compounds [1]. In a BN, the state is defined as a vector of Boolean variables, each representing a distinct component of the system under consideration. The time evolution of each variable is governed by an update function, i.e., a Boolean expression that encodes how the other variables affect the change of state at each (discrete) time step [2,3,4,5]. In the main text we focus on synchronous BNs whereby the next state is obtained by applying all the update functions to the activation values of the current state. However, in the supplementary material we show how our approach can be applied also to BNs with partially asynchronous update schema (see, e.g., the priority classes supported by the popular tool GINsim [6] as described in [7]).
From a computational viewpoint, BNs are challenging to analyze. For example, the state space of the network, known as the state transition graph (STG), has exponential size in the number of variables. Thus, a full enumeration of the state space is possible only for networks of limited size. Another relevant type of analysis concerns the computation of attractors, i.e., those sets of states toward which the system tends to evolve and remain [8, 9]; these are often associated with biologically intelligible conditions of the system under study such as cell differentiation [3, 10]. Attractor identification is NPhard [11] and, even if efficient tools have been developed [12], they do not scale well for large BNs.
These computational difficulties have motivated the development of reduction methods to ease BN analysis. Available techniques can be classified in three families according to the type of reduction: (i) by reasoning directly on the BN structure [4, 13,14,15,16]; (ii) by reducing the underlying STG [5, 17]; (iii) by transforming a BN into other formalisms for which specific reduction techniques are available [18, 19]. The latter two classes suffer two main limitations. First, STGbased reductions are still subject to state space explosion since they require the full enumeration of the state space to start with. Second, reductions via other formalisms may not be complete in the sense that the dynamics of the original BN and of the transformed model are not equivalent, hence some reductions may be missed (see Additional file 1).
In the case of reduction methods at the BN level, popular examples are based on the notion of variable absorption, proposed originally in [14, 15]. The main idea is that certain BN variables can get removed by replacing their occurrences with their update functions. This is based on the assumption that those variables evolve over time scales that justify that they can be updated first in the model. Other methods remove output/leaf variables [4, 13] (variables that do not appear in the update functions of other variables) or frozen ones (variables that stabilize to the same value after some iterations independently of the initial conditions) [16].
Here we present a complementary type of reduction method based on the computation of a partition of the variables in the BN, whereby the future dynamics of variables in a block of the partition are equal whenever they start from the same condition. This can be convenient, for instance, if one is interested in studying the dynamics due to simultaneous activation or deactivation of groups of variables [20] (see also the case studies presented in the Results and discussion section). We call this kind of relation a Boolean backward equivalence (BBE) because it is defined analogously to the notion of backward bisimulation for Markov chains [21], more recently extended to chemical reaction networks [19, 22] and ordinary differential equations [23]. Recently, it has been shown how this backward notion can be given also for linear differential algebraic equations (DAE) [24, 25]. Using DAE terminology, such backward notion has been related to the preservation of invariant subspaces.
Every reduction technique comes with its own intuitive interpretation. For example, if we consider variable absorption mentioned above, it is intuitively based on the idea of fast/slow decomposition: it is biologically plausible to absorb a variable in the update function of another when the former fires faster than the latter. BBE is based on the following three orthogonal considerations:

BBE allows the modeler to discover chains of variables that, under some initialization conditions, describe the same dynamics. This might be interesting, e.g., to estimate the quality of a model: large BBE reductions might signal excessive redundancy in the model.

As mentioned above, the modeler might be interested only in dynamics where two or more variables have simultaneous (de)activation value (see, e.g., [20]). The TLGL case study in section Results and discussion further discusses this.

In [23], it has been shown that this backward notion corresponds to Cardelli’s emulation [26] which enables to relate a complex model with a simpler one. Interestingly, [26] discusses how emulation can be given an evolutionary interpretation. In fact, an original model can express all the dynamics of the reduced model. In addition, the original model can also express all additional dynamics coming from states where variables related by emulation have different activation values (not permitted in the reduced model because variables related by emulation get collapsed in the same reduced one). Given this richer dynamics of the original model, Cardelli uses selected case studies in [26] to argue how the original model can be seen as an evolved version of the reduced one. We do not further investigate this aspect for BBE. However, given that BBE is based as well on the mentioned backward notions, it is not surprising that there exists a similar relation among the dynamics expressed by the original and reduced model (cf Fig. 1).
The criteria for a candidate partition of variables to be a BBE are encoded into a satisfiability problem over the expressions of the BN’s update functions: we synthesise a Boolean expression involving BN variables and check whether there exists at least one combination of truth values for the variables that makes such expression true. This type of test can be effectively implemented using tools known as SAT solvers [27].
If a partition is a BBE, a reduced BN can be obtained by choosing and maintaining only a representative variable for each partition block, and renaming all variables in the remaining update functions with the representative one from their block. The STG of the reduced network exactly preserves the original dynamics for all states that have equal values across variables in the same block (Fig. 1). Importantly, however, the reduction method does not require the generation of the original STG, making it possible to obtain a reduced STG also from instances that would not be analyzable due to their massive size.
A crucial property satisfied by BBE is that there exists a maximal reduction for each BN, i.e., the coarsest BBE partition. This can be computed using a partitionrefinement algorithm in a similar fashion as in Markov chains [28], reaction networks [22] and differential equations [23]. The algorithm essentially builds upon a fundamental result in computer science to prove equivalences in formal languages [29]. Given an initial partition of variables, the algorithm splits the blocks of the partition to compute its coarsest refinement that satisfies the BBE criteria. Thus, the maximal reduction is obtained when all variables are in the same unique block of the initial partition. However, the possibility of arbitrarily choosing the initial partition unlocks modelspecific reduction queries that preserve the dynamics of userdefined variables. For example, in typical BN models of signaling pathways [30, 31], certain variables may represent the input signals to upstream components such as receptors. Formally, inputs may be detected because their update functions are constants that represent the values of such inputs. In this case, a possibly more biologically relevant initial partition may separate inputs from the other variables, obtaining inputseparated (IS) reductions.
Our partitionrefinement algorithm takes a polynomial number of steps as a function of the number of BN variables. At each iteration, it queries a SAT solver to check for the BBE criteria. If the query is satisfiable, i.e., the current partition is not a BBE, the returned assignment is used to split the current partition and perform another iteration; if the query is unsatisfiable, it returns the current partition as the coarsest BBE refinement of the initial one. We fully develop our approach in the Methods section and in Additional file 2.
Interestingly, although our algorithm is theoretically as complex as SAT solving, it behaves effectively in practice. Using a prototype implementation available within the software tool ERODE [32], we demonstrate its performance on a largescale validation across 86 BN models from two wellknown repositories [6, 33]. We show that almost all BNs can be reduced by BBE, providing speedups for the computation of STGs and attractors by more than three orders of magnitude. In some cases, BBE could render the analysis feasible in instances that originally issued outofmemory errors or that were stopped after long time outs. This comes at the cost that part of the original dynamics is lost. In particular, in the STG we preserve all and only the states where variables within the same BBEblock have the same value, and transitions among them. From this, and from the properties of BBE, we also get that the method preserves all and only the attractors containing at least one preserved state. This confirms that BBE is complementary to existing reduction techniques for BNs. Indeed, in several areas of science and engineering, it is common to have reduction techniques that:

Preserve all dynamics but might add spurious ones. An example is [14] which preserves all attractors but might create new spurious ones. These often come with the name of overapproximations (this is because, e.g., [14] might overapproximate the set of attractors of a model by computing a larger set containing all original ones, plus some spurious ones;

Do not preserve all the dynamics, but guarantee to not add spurious ones, like BBE. These often come with the name of underapproximations (this is because, e.g., BBE might underapproximate the set of attractors of a model by computing a smaller set containing only original ones, but potentially not all).
These two families of techniques are not comparable. They might be jointly used to obtain upperbounds (the case of [14]), and lowerbounds (the case of BBE) on the actual number of attractors in a model.
This paper extends the previous conference version [34]. All numerical experiments have been redesigned by adding an additional model repository, by performing a largescale validation of the analysis speedups offered by BBE, and by considering a more recent and efficient tool for identification of attractors. We have also performed a new largescale validation on randomly generated BNs. Finally, we have generalised the theory to support also BNs with partially asynchronous update schema, and we have added in Additional file 3 a new case study considering one such BN.
Results and discussion
In this section, we perform a largescale validation of BBE. We first check its reduction power on published models from the literature, and then we demonstrate how it facilitates the analysis tasks of STG generation and attractor computation. In particular, we show how BBE brings important analysis speedups, both in terms of STG generation, and attractor analysis. In several cases, BBE enables the analysis of models that were originally intractable due to their complexity.
After this, we use two selected case studies to show how one can tune the reduction power of BBE to preserve or exclude specific dynamics of interest.
Toolchain. We implemented our method in ERODE [32], a freely available software for the modeling, analysis, and reduction of biological systems modelled in terms of BNs [34], differential equations [35], and chemical reaction networks [19]. ERODE integrates the SAT solver Z3 [36]. Thanks to importing/exporting functionalities, we let ERODE interact with the COLOMOTO Notebook [37, 38], which integrates several tools for modeling and analysis of BNs. STG generation is performed using the tool PyBoolNet [39], while attractor identification is performed using the SATbased tool BNS [12].^{Footnote 1}
Configuration. All experiments were conducted on a machine equipped with an Intel Xeon(R) 2.80 GHz processor and 32 GB RAM. We imposed an arbitrary timeout of 8 h for each task, after which we terminated the analysis. We refer to these cases as timeout, while we use outofmemory if the execution issued a memory error.
We conducted our investigation using two model repositories: GINsim repository [6] (http://ginsim.org/models_repository), which contains 83 models, and the Biomodels repository [42] (https://www.ebi.ac.uk/biomodels/), which contains 24 models, obtaining overall 98 distinct models (9 appeared in both repositories). From these, we restricted only to models with input variables, obtaining 86 models. In other words, we considered about 92% of the models available in the two respositories. This selection was done to avoid favouring BBE: in BNs without inputs, IS initial partitions correspond to the maximal ones, which, as the name says, allow for the best possible BBE reduction of a model in terms of aggregation power. Part of these 86 models, 45, are multivalued networks, i.e. logical models wherein some variables take more than two activation statuses, e.g., \(\{0,1,2\}\) for low, medium, or high concentration respectively (see, e.g., [43]). We transform such models in dynamically equivalent BNs by applying a socalled booleanization technique [44], supported by GINsim [40].
We consider two reduction scenarios relevant to input variables, using maximal and inputseparated (IS) initial partitions like \(\mathcal {H}_0\) in Eq. (1) and \(\mathcal {H}_0'\) in Eq. (2), respectively. In Additional file 4 we perform a similar analysis on randomly generated BNs.
Largescale validation
Largescale validation: reduction power. We begin by addressing the reduction power of BBE. For this, we consider the reduction ratios (variables in the reduced BN over the variables in the original one) obtained on all models.
Figure 2 displays the reduction ratios for both the maximal and IS reductions. We observe that almost all models can be reduced by BBE, in particular 93% admit a maximal reduction, while 91% admit an IS one. The reduction ratios distribute almost uniformly from 0.15 to 1.00 (no reduction), with average reduction ratio of 0.70 and 0.77 for maximal and IS reductions, respectively. For most models, the maximal and IS reduction ratios do not change significantly, meaning that BBE is effective also when we prevent input variables from merging with internal variables. All detailed results of this analysis can be found in Additional file 5: Table S3.
Largescale validation: STG generation speedup. We hereby demonstrate the speedups that BBE provides to STG generation on a selection of the considered models. Figure 3 focuses on the 20 models with more than 10 variables for which the STG generation succeeded in both the original and reduced models, while Fig. 4 focuses on the 13 ones where the STG generation failed on the original model and succeeded in the maximal or IS reduction. On the used machine, STG generation failed for models with 24 or more variables. Therefore, Fig. 3 focuses on models with less than 24 variables, while Fig. 4 focuses on models with 24 or more variables for which at least one reduction had less than 24 variables.
The red line in Fig. 3 marks the area where the reduction would not bring speedups. We can see that all points are below the line, with instances showing more than two orders of magnitude difference between the original and reduced runtimes. This proves that BBE can effectively lead to faster STG generation. All cases where the dots and crosses overlap refer to models where the two reductions coincide.
We now consider the 13 models in Fig. 4 where STG generation was not feasible for the original models. We note that the generation succeeded for all maximal reductions, while it failed for two IS ones. As denoted by the model identifiers in the xaxis, these are models 4 and 15 in Fig. 2, where the IS reductions have 34 and 25 variables, respectively. The largest runtime is 441 s for the IS reductions, and 338 s for the maximal ones. Detailed results are presented in Additional file 5: Table S4.
Largescale validation: Attractor computation speedup Fig. 5 studies the speedups that BBE provided to the computation of attractors on the models from Fig. 2. The plot has the same structure as Fig. 3. We observe that, in several cases, we have significant analysis speedups. In particular, we note how the dots and crosses spread to the right, due to original runtimes in the order of \(10^3\,\text {s}\), while they hardly go up beyond \(1\,\text {s}\) for runtimes on reduced models. Furthermore, models 18 and 29 from Fig. 2 are omitted here because the analysis failed on the original models. Instead, the analysis of their maximal reductions required at most 0.15 s, and that of their IS reductions at most 2.5 s. Detailed results are given in Additional file 5: Table S5.
Largescale validation: Interpretation. BBE can successfully reduce a large amount of models. For the original models, the state space explosion prevents full state space exploration in many cases, and hampers the identification of the attractors. This is mitigated in practice by BBE, with extreme cases where BBE made analysis feasible whereas the original models were intractable. As shown in Additional file 5: Table S5, part of the attractors are lost in the reduced models, namely those not involving constant states on the computed BBE (see Methods section). Table S5 shows cases like model 18 or 29, whose attractors could not be computed at all without BBE reduction. At the same time, the table shows that, by using the default IS or maximal initial partitions, a large part of the attractors might be lost (because, according to our theory, involve STG states where variables belonging to the same BBE block have different value). In Additional file 4: Fig. S7 we provide this information graphically for the BNs from the two repositories, and for randomly generated ones. This can be mitigated by devising refined initial partitions for the model and problem at hand. This is exemplified and discussed in greater detail in the next section where we show how a modeler can easily devise refined initial partitions that may allow to preserve more attractors, or drop those that are not of interest.
Case studies
In the previous part of this section we studied the aggregation power and the analysis speedups offered by BBE on 86 models from the literature. Here, instead, we use two selected case studies (MAPK, TLGL) to show how one can tune, or refine the reduction power of BBE using modelspecific information to preserve all dynamics of interest, and selectively exclude behavior that does not have biological relevance. Nevertheless, for completeness, we provide in Table 1 information on the analysis runtimes on all models (and their reductions) discussed in this section. We consider analysis runtimes on the models (Original), and on their IS (IS) and maximal (Maximal) reductions as done in the largescale validation. Furthermore, we consider an additional reduction obtained using a refined initial partition discussed in the corresponding sections (Refined). STG generation failed on all models and reductions because they all have more than 24 variables. Indeed, we have previously discussed how, on the used machine, STG generation fails for models with 24 or more variables. Instead, attractor analysis succeeded on all models, with important speedups obtained for all reductions. For both models, the IS and Maximal cases have a particularly low analysis runtime. This is because, as we shall discuss next, several attractors are discarded in these reductions. Notably, despite the Refined reductions have speedup factors of about two, as we shall see they preserve all attractors for MAPK, and all attractors of interest for TLGL.
MAPK case study We consider a BN model for MitogenActivated Protein Kinase (MAPK) from [45]. The model consists of tightly interconnected signalling pathways involved in diverse cellular processes, such as cell cycle, survival, apoptosis and differentiation. The BN is depicted in Fig. 6. It contains 53 variables, 4 of which being inputs (\(EGFR\_stimulus\), \(FGFR3\_stimulus\), \(TGFBR\_stimulus\), and \(DNA\_damage\)), and has 40 attractors.
MAPK: Maximal and IS reduction The maximal BBE reduction of this model has 39 variables. The discovered blocks are visualized in Fig. 6 using different background colors. In particular, we note that the yellow block contains all inputs and five noninputs variables, three related to \(TGFBR\_stimulus\), and two related to \(DNA\_damage\). Instead, the IS reduction has 41 variables, the only difference being that the block with inputs from the maximal reduction (Fig. 6) is split in three blocks: one for the inputs, one for the two noninput variables directly connected to the two rightmost inputs, and one for the remaining noninput variables. In both cases, the reduced BNs have 17 attractors.
MAPK: Refined reduction with vertical IS We propose a third modelspecific initial partition that considers inputs also indirectly. Intuitively, variables like \(TGFBR\) depend only on the value assigned to an input (\(TGFBR\_stimulus\)). This reasoning can be iterated downward through the pathway, allowing to add also \(TAK1\), and \(SMAD\), until variables that depend on other (input) variables are met. In some sense, we can see \(TGFBR\), \(TAK1\), and \(SMAD\) as indirect inputs. This is because, in a few iterations the value assigned to the corresponding input will be propagated to them, and they will not change value anymore. In other words, we use a block per input, each containing the input and all noninput variables only positively affected by the input or by variables in the block. That way, we obtain an initial partition denoted by the blue dashed shapes in Fig. 6, plus an additional fifth block containing all other variables. The rationale is that a variable only affected by an input will have the same truth value of the input, therefore it can be considered as a sort of indirect input. The obtained BBE is depicted in Fig. 7. The reduced BN contains 42 variables and preserves all 40 attractors.
TLGL case study We consider a BN model for TLGL from [20]. It refers to the disease TLGL leukemia which features a clonal expansion of antigenprimed, competent, cytotoxic T lymphocytes (CTL). This BN is a signalling pathway, constructed empirically through extensive literature review, and determines the survival of CTL. The BN, depicted in Fig. 8, consists of 60 variables, 6 of which are inputs (the yellow nodes in Fig. 8). The model has 264 attractors.
TLGL: Maximal and IS reduction The maximal and IS BBE coincide, as depicted in Fig. 8. We have only two nonsingleton blocks: one consisting of all the inputs, and one consisting of \(FasT\), \(A20\), \(TNF\), and \(RANTES\). The reduced BN has 52 variables and 6 attractors, which means that most of the attractors are lost.
TLGL: Refined reduction In [20], the authors discover that the simultaneous activation of the two input variables \(IL15\) and \(PDGF\) is sufficient to produce all dynamics of interest to them (namely, all the known socalled deregulations and signalling abnormalities).
In terms of initial partitions for BBE, we can encode the notion of contemporary activation or deactivation of the two inputs by using a model and problemspecific initial partition where \(IL15\) and \(PDGF\) form a block. Furthermore, we assign every input to a singleton block, while all noninput variables belong to the same block, for a total of 56 blocks. It turns out that this initial partition is actually a BBE, which therefore does not get refined by our algorithm. The reduced BN has 120 attractors.
Conclusion
Boolean backward equivalence (BBE) is an automatic reduction technique for Boolean networks (BNs) which exactly preserves dynamics of interest to the modeler by collapsing variables that are proven to have equal values in all states. The method, based on a partition refinement algorithm, can be tuned on a model and problemspecific way by specifying which variables should be preserved using an appropriate choice of the initial partition. The approach is complementary to the state of the art. Roughly, in [4, 13], reduction is achieved by replacing variables with constants and propagating those in the transitions of somehow richer STGs or across the network, respectively. Thus, the reduced model cannot be used to investigate how changes in those variables affect the dynamics. In a BBE reduction, instead, variables are collapsed into blocks and the original dynamics is exactly recovered whenever variables in the same block are assigned equal values. These studies [4, 13] additionally remove the output [4] variables (also called leaf variables [13]). However, output variables sometimes are used to denote different “responses” by the modelled system [4, 30], therefore their removal might not always be appropriate.
In variable absorption [14, 15], the main assumption is that there are variables that are updated faster than others, therefore one class of variables can be assumed to be constant and absorbed if focusing on the dynamics of the other class. Unlike BBE, this can only increase the number of attractors. In particular, variable absorption preserves exactly all steady states (singlestate attractors), while it might change the length of other attractors. Furthermore, new spurious attractors might be added. Instead, BBE might decrease the number of attractors (it discards all and only the attractors involving states where BBEequivalent variables have different activation values), but all preserved attractors are preserved exactly, including their length and reachibility from (preserved) initial states, and no spurious ones are added. Regarding other relevant work, in [16], the authors identify variables that have the same value in attractors only, but, differently from BBE, might behave differently in other states of the STG.
We validated BBE on 86 BNs from two model repositories, providing reductions and analysis speedups in almost all cases. In some, BBE enabled the analysis of models which would be otherwise intractable. There were also instances for which the reduced model could not be analyzed. This calls for further research into more aggressive reductions; for example, in its current implementation multivalued BNs are first translated into ordinary BNs, but this causes a blowup in the number of variables. It is worth investigating approaches that circumvent the intermediate translation to reduce dimensionality. Another area of research concerns the different semantic interpretations of a BN. Currently, BBE supports BN with synchronous and partially asynchronous updates; we plan to investigate variants of BBE for probabilistic BNs.
Methods
Here we explain the key steps of the reduction procedure on the BN in Fig. 9. Its customary graphical representation allows one to distinguish different kinds of variables depending on whether they appear in the update functions of other variables (as indicated by the green arrows). In the example, TLR5 can be interpreted as an input because its state remains constant and unaffected by other variables. Inputs, which often denote external stimuli [4], are explicitly set by the modeler to perform experiment campaigns. Conversely, IRAK4 and PIK3AP1 can be considered output variables because they do not appear in the update functions of other variables.
Step 1: initial partition. Our reduction algorithm starts with the specification of an initial partition of variables. The idea behind initial partitions is that the modeler can force our algorithm to not collapse given variables, by placing them in different initial blocks. In the case studies presented in the Results and discussion section we see examples of userspecified initial partitions enabling analyses of interest on the considered models. This is how initial partitions shall be used, devising casebycase useful ones. In order to favour a systematic largescale validation of our approach, here we consider two examples of initial partitions whose computation can be easily automated: the maximal partition, where all variables are placed in the same block; and the inputseparated (IS) one, where the inputs are separated from the other variables, i.e., we use an initial partition with two blocks, one for input variables and one for the other variables. ,^{Footnote 2} In the example, these are respectively given by the partitions:
and
Iterative step: splitting by the BBE condition. At every iteration, the algorithm checks the BBE condition on the current partition. Formally, BBE is defined as a partition \(X\) of variables that renders the following formula valid:
This is a Boolean formula for: whenever all variables in the same block have same value, they will not be distinguished in the next state. In other words, \(\Phi ^{\mathcal {H}}\) says that if for all partition blocks \(H_i\) the variables in \(H_i\) are equal, then the evaluations of update functions of variables in the same block stay equal. A SAT solver can determine if \(\Phi ^{\mathcal {H}}\) is valid by checking the unsatisfiability of its negation. For example, given the \(\mathcal {H}'_0\) partition in Eq. 2, one can obtain that \(\lnot \Phi ^{\mathcal {H}'_0}\) is satisfiable (i.e., \(\mathcal {H}'_0\) is not a BBE) because there exists the assignment s given by
for which, as it can be seen in the STG of Fig. 9, the next state \(s'\) is
This assignment proves that variables \(x_{ TICAM1 }\), \(x_{ MyD88 }\), \(x_{ IRAK4 }\), and \(x_{ PIK3AP1 }\) cannot belong to the same block of a partition that satisfies the BBE criteria because despite having the same value (0) in the source state s, they differ in the target state \(s'\). In addition, the assignment s suggests to split that block into two subblocks for which that assignment does not disprove the BBE condition: \(x_{ TICAM1 }\) and \(x_{ MyD88 }\) have same value in \(s'\), as well as \(x_{ IRAK4 }\) and \(x_{ PIK3AP1 }\). Thus the algorithm will perform a new iteration with the refined partition
With this, \(\lnot \Phi ^{\mathcal {H}'_1}\) is unsatisfiable, implying that \(\mathcal {H}'_1\) is a BBE partition. In Theorem 2 from Additional file 2, we prove that this algorithm returns, for any initial partition, its unique coarsest refinement that satisfies the BBE condition (3). Overall, the algorithm takes at most n steps, where n is the number of BN variables; at every step, it iterates through the provided SAT assignment, if any is provided, to perform the splitting. Thus, overall the algorithm is as hard as SAT solving; however, the numerical evaluation presented in the Results and discussion section will show how it can effectively tackle BN models from the literature.
BBE properties As discussed in Fig. 1, given a BBE it is possible to construct a reduced BN where each variable represents a partition block (Proposition 4 from Additional file 2). The STG of the reduced BN agrees with the original STG on all, and only on, states that are constant on the partition, i.e., whose variables in the same block have the same value (Proposition 4 from Additional file 2). The reduction also preserves any attractor of the original BN which contains at least one state that is constant on the partition (Theorem 5 from Additional file 2). Thus, in particular the reduced BN maintains the exact length of the attractors that are preserved without introducing spurious dynamical behavior. Instead, all states non constant on the partition are dropped, as well as all attractors not containing any state constant on the BBE partition.
We use two examples to better explain the exact preservation of part of the attractors. Considering preserved attractors, we have seen in Fig. 1 that the twostates attractor of the original model (Fig. 1 bottomleft) is preserved in a twostates attractor in the reduced model (Fig. 1 bottomright). This is the case for any preserved attractor; the number of states is preserved.
As regards attractors that are not preserved, we provide in Fig. 10 (topleft) a simple BN with 3 variables (\(x_1\), \(x_2\), and \(x_3\)) and 4 attractors (steadystates, Fig. 10, bottomleft). Figure 10 (topright) shows a BBE reduction of the model where \(x_1\) and \(x_2\) get collapsed. We can see in Fig. 10 (bottomright) that 2 attractors are preserved in the BBE reduction, while the other 2 attractors belong to the part of the STG that is not preserved, and therefore are not present in the reduced BN. In particular, according to our theory, the two attractors where \(x_1\) and \(x_2\) are both 1 or both 0 are preserved. Instead, the other two attractors have different values for \(x_1\) and \(x_2\), and therefore are not preserved.
Partially asynchronous BNs In Additional file 2, we show how BBE can also be applied to partially asynchronous BNs. Here, we equip a BN with a partition \(\mathcal {K}\) of its variables that we name synchronization partition. A new state is obtained by selecting one of the blocks K of \(\mathcal {K}\), and then applying of the update functions of the variables in K only. The activation values of the other variables are not modified. Notably, this synchronization schema is supported, e.g., by popular BN analysis tools like GINsim [6] under the name of priority classes [7].^{Footnote 3} In particular, BBE can be applied to such BNs with the caveat that the initial partition must be \(\mathcal {K}\) or refinements of it. In Additional file 3 we apply BBE to a BN with partially asynchronous update schema.
Availability of data and materials
Material to replicate the largescale experiments, including all models, are available here: https://www.erode.eu/models/BMCBioInf_CMSB2021.zip
Notes
Among the tools available in the COLOMOTO notebook, we could have opted for GINsim [40] and BoolSim [41] for STG generation and attractor identification, respectively, for BNs with synchronous update schema. In both cases, we have opted for the tools with best performances according to preliminary experiments we conducted, not reported here.
We refer to [34] for a third example of initial partition, inputdistinguished where inputs were further separated from each other. As exemplified in the case studies in the Results and discussion section, initial partitions should be defined by the modeler depending on the model at hand and on the properties to be studied. We discuss the maximal and IS partitions here to enable the largescale validation of BBE discussed in the same section.
The dynamics of BNs considered in this paper are less general than those offered by GINsim using priority classes, as they also further allow to assign different priorities to the classes, and to update the variables within them asynchronously.
Abbreviations
 BN:

Boolean network
 BBE:

Boolean backward equivalence
 STG:

State transition graph
References
Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969;22(3):437–67.
Wang RS, Saadatpour A, Albert R. Boolean modeling in systems biology: an overview of methodology and applications. Phys Biol. 2012;9(5):055001.
Azpeitia E, Benítez M, Vega I, Villarreal C, AlvarezBuylla ER. Singlecell and coupled GRN models of cell patterning in the Arabidopsis thaliana root stem cell niche. BMC Syst Biol. 2010;4(1):1–19.
Naldi A, Monteiro PT, Chaouiya C. Efficient handling of large signallingregulatory networks by focusing on their core control. In: International conference on computational methods in systems biology. Springer; 2012. p. 288–306.
Bérenguier D, Chaouiya C, Monteiro PT, Naldi A, Remy E, Thieffry D, Tichit L. Dynamical modeling and analysis of large cellular regulatory networks. Chaos Interdiscip J Nonlinear Sci. 2013;23(2):025114.
Naldi A, Berenguier D, Fauré A, Lopez F, Thieffry D, Chaouiya C. Logical modelling of regulatory networks with GINsim 2.3. Biosystems. 2009;97(2):134–9.
Fauré A, Naldi A, Chaouiya C, Thieffry D. Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle. Bioinformatics. 2006;22(14):124–31.
Schwab JD, Kühlwein SD, Ikonomi N, Kühl M, Kestler HA. Concepts in Boolean network modeling: What do they all mean? Comput Struct Biotechnol J. 2020;18:571–82. https://doi.org/10.1016/j.csbj.2020.03.001.
Hopfensitz M, Müssel C, Maucher M, Kestler HA. Attractors in Boolean networks: a tutorial. Comput Stat. 2013;28(1):19–36.
Drossel B. Random Boolean networks. Rev Nonlinear Dyn Complex. 2008;1:69–110.
Akutsu T, Kuhara S, Maruyama O, Miyano S. A system for identifying genetic networks from gene expression patterns produced by gene disruptions and overexpressions. Genome Inform. 1998;9:151–60.
Dubrova E, Teslenko M. A satbased algorithm for finding attractors in synchronous Boolean networks. IEEE/ACM Trans Comput Biol Bioinform. 2011;8(5):1393–9.
Bilke S, Sjunnesson F. Stability of the Kauffman model. Phys Rev E. 2001;65(1): 016129.
Naldi A, Remy E, Thieffry D, Chaouiya C. Dynamically consistent reduction of logical regulatory graphs. Theor Comput Sci. 2011;412(21):2207–18.
VelizCuba A. Reduction of Boolean network models. J Theor Biol. 2011;289:167–72.
Richardson KA. Simplifying Boolean networks. Adv Complex Syst. 2005;8(04):365–81.
Figueiredo D. Relating bisimulations with attractors in Boolean network models. In: International conference on algorithms for computational biology. Springer; 2016. p. 17–25.
Zañudo JGT, Albert R. An effective network reduction approach to find the dynamical repertoire of discrete dynamic networks. Chaos Interdiscip J Nonlinear Sci. 2013;23(2): 025111. https://doi.org/10.1063/1.4809777.
Cardelli L, Tribastone M, Tschaikowski M, Vandin A. Maximal aggregation of polynomial dynamical systems. Proc Natl Acad Sci. 2017;114(38):10029–34.
Zhang R, Shah MV, Yang J, Nyland SB, Liu X, Yun JK, Albert R, Loughran TP. Network model of survival signaling in large granular lymphocyte leukemia. Proc Natl Acad Sci. 2008;105(42):16308–13.
Sproston J, Donatelli S. Backward bisimulation in Markov chain model checking. Softw Eng IEEE Trans. 2006;32(8):531–46. https://doi.org/10.1109/TSE.2006.74.
Cardelli L, Tribastone M, Tschaikowski M, Vandin A. Forward and backward bisimulations for chemical reaction networks. In: 26th International conference on concurrency theory, CONCUR 2015, Madrid, Spain, September 1.4, 2015. 2015. p. 226–239. https://doi.org/10.4230/LIPIcs.CONCUR.2015.226.
Cardelli L, Tribastone M, Tschaikowski M, Vandin A. Symbolic computation of differential equivalences. ACM SIGPLAN Not. 2016;51(1):137–50.
Tognazzi S, Tribastone M, Tschaikowski M, Vandin A. Differential equivalence for linear differential algebraic equations. IEEE Trans Autom Control. 2022;67(7):3484–93. https://doi.org/10.1109/TAC.2021.3108530.
Tognazzi S, Tribastone M, Tschaikowski M, Vandin A. Backward invariance for linear differential algebraic equations. In: 2018 IEEE conference on decision and control (CDC). 2018. p. 3771–3776. https://doi.org/10.1109/CDC.2018.8619710
Cardelli L. Morphisms of reaction networks that couple structure to function. BMC Syst Biol. 2014;8(1):84.
Biere A, Biere A, Heule M, van Maaren H, Walsh T. Handbook of satisfiability: volume 185 frontiers in artificial intelligence and applications. IOS Press, NLD. 2009.
Valmari A, Franceschinis G. Simple O(m logn) time markov chain lumping. In: Tools and algorithms for the construction and analysis of systems, 16th international conference, TACAS 2010, held as part of the joint European conferences on theory and practice of software, ETAPS 2010, Paphos, Cyprus, March 2028, 2010. Proceedings. 2010. p. 38–52. https://doi.org/10.1007/9783642120022_4.
Paige R, Tarjan R. Three partition refinement algorithms. SIAM J Comput. 1987;16(6):973–89. https://doi.org/10.1137/0216062.
RodríguezJorge O, KempisCalanis LA, AbouJaoudé W, GutiérrezReyna DY, Hernandez C, RamirezPliego O, ThomasChollier M, Spicuglia S, Santana MA, Thieffry D. Cooperation between T cell receptor and tolllike receptor 5 signaling for CD4+ T cell activation. Sci Signal. 2019. https://doi.org/10.1126/scisignal.aar3641.
Klamt S, SaezRodriguez J, Lindquist JA, Simeoni L, Gilles ED. A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinform. 2006;7(1):56.
Cardelli L, Tribastone M, Tschaikowski M, Vandin A. Erode: a tool for the evaluation and reduction of ordinary differential equations. In: International conference on tools and algorithms for the construction and analysis of systems. Springer; 2017. p. 310–328.
Le Novere N, Bornstein B, Broicher A, Courtot M, Donizelli M, Dharuri H, Li L, Sauro H, Schilstra M, Shapiro B, et al. Biomodels database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems. Nucleic Acids Res. 2006;34(suppl–1):689–91.
Argyris G, Lluch Lafuente A, Tribastone M, Tschaikowski M, Vandin A. Reducing boolean networks with backward Boolean equivalence. In: International conference on computational methods in systems biology. Springer; 2021. p. 1–18.
Cardelli L, Tribastone M, Tschaikowski M, Vandin A. Symbolic computation of differential equivalences. In: POPL 2016. 2016. p. 137–150. https://doi.org/10.1145/2837614.2837649
De Moura L, Bjørner N. Z3: an efficient SMT solver. In: International conference on tools and algorithms for the construction and analysis of systems. Springer; 2008. p. 337–340.
Naldi A, Hernandez C, Levy N, Stoll G, Monteiro PT, Chaouiya C, Helikar T, Zinovyev A, Calzone L, CohenBoulakia S, Thieffry D, Paulevé L. The CoLoMoTo interactive notebook: accessible and reproducible computational analyses for qualitative biological networks. Front Physiol. 2018;9:680. https://doi.org/10.3389/fphys.2018.00680.
Naldi A, Monteiro PT, Müssel C, Kestler HA, Thieffry D, Xenarios I, SaezRodriguez J, Helikar T, Chaouiya C, for Logical Models C, Tools. Cooperative development of logical modelling standards and tools with CoLoMoTo. Bioinformatics. 2015;31(7):1154–9.
Klarner H, Streck A, Siebert H. PyBoolNet: a python package for the generation, analysis and visualization of Boolean networks. Bioinformatics. 2017;33(5):770–2.
Chaouiya C, Naldi A, Thieffry D. Logical modelling of gene regulatory networks with GINsim. 2012;804:463–79.
Di Cara A, Garg A, De Micheli G, Xenarios I, Mendoza L. Dynamic simulation of regulatory networks using squad. BMC Bioinform. 2007;8(1):462.
...MalikSheriff RS, Glont M, Nguyen TVN, Tiwari K, Roberts MG, Xavier A, Vu MT, Men J, Maire M, Kananathan S, Fairbanks EL, Meyer JP, Arankalle C, Varusai TM, KnightSchrijver V, Li L, DueñasRoca C, Dass G, Keating SM, Park YM, Buso N, Rodriguez N, Hucka M, Hermjakob H. BioModels—15 years of sharing computational models in life science. Nucleic Acids Res. 2020;48(D1):407–15. https://doi.org/10.1093/nar/gkz1055.gkz1055.
Fauré A, Vreede B, Sucena E, Chaouiya C. A discrete model of Drosophila eggshell patterning reveals cellautonomous and Juxtacrine effects. PLoS Comput Biol. 2014;10:1003527. https://doi.org/10.1371/journal.pcbi.1003527.
Delaplace F, Ivanov S. Bisimilar Booleanization of multivalued networks. BioSystems, 2020;104205
Grieco L, Calzone L, BernardPierrot I, Radvanyi F, KahnPerles B, Thieffry D. Integrative modelling of the influence of MAPK network on cancer cell fate decision. PLoS Comput Biol. 2013;9(10):1003286.
Coolen M, Thieffry D, Drivenes Ø, Becker TS, BallyCuif L. miR9 controls the timing of neurogenesis through the direct inhibition of antagonistic factors. Dev Cell. 2012;22(5):1052–64.
Wittmann DM, Krumsiek J, SaezRodriguez J, Lauffenburger DA, Klamt S, Theis FJ. Transforming Boolean models to continuous models: methodology and application to Tcell receptor signaling. BMC Syst Biol. 2009;3(1):98. https://doi.org/10.1186/17520509398.
Sipser M. Introduction to the theory of computation. 3rd ed. Boston: Course Technology; 2013.
Müssel C, Hopfensitz M, Kestler HA. BoolNet–an R package for generation, reconstruction and analysis of Boolean networks. Bioinformatics. 2010;26(10):1378–80.
Mbodj A, Junion G, Brun C, Furlong EE, Thieffry D. Logical modelling of Drosophila signalling pathways. Mol BioSyst. 2013;9(9):2248–58.
MartinezSanchez ME, Hiriart M, AlvarezBuylla ER. The CD4+T cell regulatory network mediates inflammatory responses during acute hyperinsulinemia: a simulation study. BMC Syst Biol. 2017;11(1):1–12.
Mbodj A, Gustafson EH, Ciglar L, Junion G, Gonzalez A, Girardot C, Perrin L, Furlong EE, Thieffry D. Qualitative dynamical modelling can formally explain mesoderm specification and predict novel developmental phenotypes. PLoS Comput Biol. 2016;12(9):1005073.
MartinezSanchez ME, Mendoza L, Villarreal C, AlvarezBuylla ER. A minimal regulatory network of extrinsic and intrinsic factors recovers observed patterns of CD4+ T cell differentiation and plasticity. PLoS Comput Biol. 2015;11(6):1004324.
Fauré A, Vreede BM, Sucena É, Chaouiya C. A discrete model of Drosophila eggshell patterning reveals cellautonomous and Juxtacrine effects. PLoS Comput Biol. 2014;10(3):1003527.
Sánchez L, Chaouiya C. Primary sex determination of placental mammals: a modelling study uncovers dynamical developmental constraints in the formation of Sertoli and granulosa cells. BMC Syst Biol. 2016;10(1):1–11.
Mombach JC, Bugs CA, Chaouiya C. Modelling the onset of senescence at the G1/S cell cycle checkpoint. BMC Genom. 2014;15(S7):7.
CorralJara KF, Chauvin C, AbouJaoudé W, Grandclaudon M, Naldi A, Soumelis V, Thieffry D. Interplay between smad2 and stat5a is a critical determinant of IL17A/IL17F differential expression. Mol Biomed. 2021;2(1):1–16.
AbouJaoudé W, Monteiro PT, Naldi A, Grandclaudon M, Soumelis V, Chaouiya C, Thieffry D. Model checking to assess Thelper cell plasticity. Front Bioeng Biotechnol. 2015;2:86.
Kondratova M, Barillot E, Zinovyev A, Calzone L. Modelling of immune checkpoint network explains synergistic effects of combined immune checkpoint inhibitor therapy and the impact of cytokines in patient response. Cancers. 2020;12(12):3600.
Vaga S, BernardoFaura M, Cokelaer T, Maiolica A, Barnes CA, Gillet LC, Hegemann B, van Drogen F, Sharifian H, Klipp E, et al. Phosphoproteomic analyses reveal novel crossmodulation mechanisms between two signaling pathways in yeast. Mol Syst Biol. 2014;10(12):767.
Naldi A, Carneiro J, Chaouiya C, Thieffry D. Diversity and plasticity of TH cell types predicted from regulatory network modelling. PLoS Comput Biol. 2010;6(9):1000912.
NuñezReza KJ, Naldi A, SánchezJiménez A, LeonApodaca AV, Santana MA, ThomasChollier M, Thieffry D, MedinaRivera A. Logical modelling of in vitro differentiation of human monocytes into dendritic cells unravels novel transcriptional regulatory interactions. Interface focus. 2021;11(4):20200061.
Terfve C, Cokelaer T, Henriques D, MacNamara A, Goncalves E, Morris MK, Iersel MV, Lauffenburger DA, SaezRodriguez J. CellNOptR: a flexible toolkit to train protein signaling networks to data using multiple logic formalisms. BMC Syst Biol. 2012;6(1):1–14.
Floc’Hlay S, Molina MD, Hernandez C, Haillot E, ThomasChollier M, Lepage T, Thieffry D. Deciphering and modelling the TGFβ signalling interplays specifying the dorsalventral axis of the sea urchin embryo. Development. 2021;148(2): 189944.
Hernandez C, ThomasChollier M, Naldi A, Thieffry D. Computational verification of large logical modelsapplication to the prediction of t cell response to checkpoint inhibitors. bioRxiv. 2020.
Fauré A, Naldi A, Lopez F, Chaouiya C, Ciliberto A, Thieffry D. Modular logical modelling of the budding yeast cell cycle. Mol BioSyst. 2009;5(12):1787–96.
Calzone L, Tournier L, Fourquet S, Thieffry D, Zhivotovsky B, Barillot E, Zinovyev A. Mathematical modelling of cellfate decision in response to death receptor engagement. PLoS Comput Biol. 2010;6(3):1000702.
Niarakis A, Bounab Y, Grieco L, Roncagalli R, Hesse AM, Garin J, Malissen B, Daëron M, Thieffry D. Computational modeling of the main signaling pathways involved in mast cell activation. Fc Recept. 2014;69–93
Sahin Ö, Fröhlich H, Löbke C, Korf U, Burmester S, Majety M, Mattern J, Schupp I, Chaouiya C, Thieffry D, et al. Modeling ERBB receptorregulated G1/S transition to find novel targets for de novo trastuzumab resistance. BMC Syst Biol. 2009;3(1):1.
MacNamara A, Terfve C, Henriques D, Bernabé BP, SaezRodriguez J. Statetime spectrum of signal transduction logic models. Phys Biol. 2012;9(4): 045003.
Selvaggio G, Canato S, Pawar A, Monteiro PT, Guerreiro PS, Brás MM, Janody F, Chaouiya C. Hybrid epithelialmesenchymal phenotypes are controlled by microenvironmental factors. Can Res. 2020;80(11):2407–20.
Chaouiya C, Bérenguier D, Keating SM, Naldi A, Van Iersel MP, Rodriguez N, Dräger A, Büchel F, Cokelaer T, Kowal B, et al. SBML qualitative models: a model representation format and infrastructure to foster interactions between qualitative modelling formalisms and tools. BMC Syst Biol. 2013;7(1):1–15.
Remy E, Rebouissou S, Chaouiya C, Zinovyev A, Radvanyi F, Calzone L. A modeling approach to explain mutually exclusive and cooccurring genetic alterations in bladder tumorigenesis. Can Res. 2015;75(19):4042–52.
González A, Chaouiya C, Thieffry D. Dynamical analysis of the regulatory network defining the dorsalventral boundary of the drosophila wing imaginal disc. Genetics. 2006;174(3):1625–34.
Sánchez L, Chaouiya C. Logical modelling uncovers developmental constraints for primary sex determination of chicken gonads. J R Soc Interface. 2018;15(142):20180165.
Hamey FK, Nestorowa S, Kinston SJ, Kent DG, Wilson NK, Göttgens B. Reconstructing blood stem cell regulatory network models from singlecell molecular profiles. Proc Natl Acad Sci. 2017;114(23):5822–9.
Béal J, Pantolini L, Noël V, Barillot E, Calzone L. Personalized logical models to investigate cancer response to BRAF treatments in melanomas and colorectal cancers. PLoS Comput Biol. 2021;17(1):1007900.
Simao E, Remy E, Thieffry D, Chaouiya C. Qualitative modelling of regulated metabolic pathways: application to the tryptophan biosynthesis in E. coli. Bioinformatics. 2005;21(suppl2):190–6.
Enciso J, Mayani H, Mendoza L, Pelayo R. Modeling the proinflammatory tumor microenvironment in acute lymphoblastic leukemia predicts a breakdown of hematopoieticmesenchymal communication networks. Front Physiol. 2016;7:349.
Sánchez L, Thieffry D. A logical analysis of the drosophila gapgene system. J Theor Biol. 2001;211(2):115–41.
Fauré A, Thieffry D. Logical modelling of cell cycle control in eukaryotes: a comparative study. Mol BioSyst. 2009;5(12):1569–81.
Mendoza L. A network model for the control of the differentiation process in TH cells. Biosystems. 2006;84(2):101–14.
González A, Chaouiya C, Thieffry D. Logical modelling of the role of the HH pathway in the patterning of the drosophila wing disc. Bioinformatics. 2008;24(16):234–40.
Sánchez L, Chaouiya C, Thieffry D. Segmenting the fly embryo: logical analysis of the role of the segment polarity crossregulatory module. Int J Dev Biol. 2002;52(8):1059–75.
Montagud A, Béal J, Tobalina L, Traynard P, Subramanian V, Szalai B, Alföldi R, Puskás L, Valencia A, Barillot E, SaezRodriguez J, Calzone L. Patientspecific Boolean models of signaling networks guide personalized treatments. bioRxiv. 2021. https://doi.org/10.1101/2021.07.28.454126.
SánchezVillanueva JA, RodríguezJorge O, RamírezPliego O, Rosas Salgado G, AbouJaoudé W, Hernandez C, Naldi A, Thieffry D, Santana MA. Contribution of ROS and metabolic status to neonatal and adult CD8+ T cell activation. PLoS ONE. 2019;14(12):0226388.
Verlingue L, Dugourd A, Stoll G, Barillot E, Calzone L, LondoñoVallejo A. A comprehensive approach to the molecular determinants of lifespan using a Boolean model of geroconversion. Aging Cell. 2016;15(6):1018–26.
Flobak Å, Baudot A, Remy E, Thommesen L, Thieffry D, Kuiper M, Lægreid A. Discovery of drug synergies in gastric cancer cells predicted by logical modeling. PLoS Comput Biol. 2015;11(8):1004426.
Zañudo JG, Steinway SN, Albert R. Discrete dynamic network modeling of oncogenic signaling: mechanistic insights for personalized treatment of cancer. Curr Opin in Syst Biol. 2018;9:1–10.
Cohen DP, Martignetti L, Robine S, Barillot E, Zinovyev A, Calzone L. Mathematical modelling of molecular pathways enabling tumour cell invasion and migration. PLoS Comput Biol. 2015;11(11):1004571.
Cacace E, Collombet S, Thieffry D. Logical modeling of cell fate specification—application to t cell commitment 2020;139:205–238
Collombet S, van Oevelen C, Ortega JLS, AbouJaoudé W, Di Stefano B, ThomasChollier M, Graf T, Thieffry D. Logical modeling of lymphoid and myeloid cell specification and transdifferentiation. Proc Natl Acad Sci. 2017;114(23):5792–9.
Traynard P, Fauré A, Fages F, Thieffry D. Logical model specification aided by modelchecking techniques: application to the mammalian cell cycle regulation. Bioinformatics. 2016;32(17):772–80.
AbouJaoudé W, Ouattara DA, Kaufman M. From structure to dynamics: frequency tuning in the P53MDM2 network: I. logical approach. J Theor Biol. 2009;258(4):561–77.
Acknowledgements
We thank reviewers of the original CMSB 2021 submission and of this submission for their help in improving the paper. We also thank the CMSB 2021 attendants for their comments and suggestions on the work. We thank Laure BallyCuif, author of [46], for her fruitful information on the modelling approach helping us in the case study in Additional file 3.
Funding
Partially supported by the DFF project REDUCTO 904000224B, the Poul Due Jensen Grant 883901, the Villum Investigator Grant S4OS, and the PRIN project SEDUCE 2017TWRCNB. The funders had no direct role in the design of the study, collection, analysis and interpretation of the data, or writing of the manuscript. This publication was produced with the cofunding European Union  Next Generation EU, in the context of The National Recovery and Resilience Plan, Investment 1.5 Ecosystems of Innovation, Project Tuscany Health Ecosystem (THE), CUP: B83C22003920001.
Author information
Authors and Affiliations
Contributions
AL, MiT, MaT, and AV conceived the study. GA, AL, MiT, MaT, and AV designed, developed and optimised the methods. GA, AL, MiT, MaT, and AV planned the experiments. GA and AV conducted the experiments and benchmarks. GA, AL, MiT, MaT, and AV wrote and revised the manuscript. All authors read and approved the final manuscript.
About this supplement
This article has been published as part of BMC Bioinformatics Volume 24 Supplement 1, 2023: Special Issue of the 19th International Conference on Computational Methods in Systems Biology. The full contents of the supplement are available online at https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume24supplement1.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1
: Comparison with encodingbased reductions.
Additional file 2
: Technical Results.
Additional file 3
: Application of BBE to a BN with partially asynchronous schema.
Additional file 4
: An application of BBE to randomly generated Boolean Networks.
Additional file 5
: Tables with detailed results from largescale experiments.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Argyris, G.A., Lluch Lafuente, A., Tribastone, M. et al. Reducing Boolean networks with backward equivalence. BMC Bioinformatics 24 (Suppl 1), 212 (2023). https://doi.org/10.1186/s12859023053269
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12859023053269
Keywords
 Boolean network
 Model reduction
 Statespace generation
 Attractors analysis
 Partition refinement