In this subsection we introduce formally the concept of percolation, used in many approaches to control. We also deduce properties of percolated subspaces that are useful for control strategy identification and that are used later in our approach.

The concept of value percolation is based on the idea that fixing of some components might induce other components to get fixed downstream. We define the percolation function with respect to f as the function that associates to every subspace S the subspace determined by the components fixed by f within S.

Definition 3.2. Given a Boolean function f, the percolation function with respect to f is the function F(f):S→S, S ↦ F(f)(S), where S is the set of all subspaces in 𝔹ⁿ and F(f)(S) is the smallest subspace that contains f(S) (with respect to inclusion).

Explicitly, given a subspace S∈S, F(f)(S) = Σ(I, c) with I = {i ∈ V | |f_i(S)| = 1} and c any state 𝔹ⁿ such that in c_i = f_i(x) for all x ∈ S for i ∈ I.

Definition 3.3. Let f be a Boolean function and S, S′ ⊆ 𝔹ⁿ two subspaces. We say that the subspace S percolates to S′ under f if and only if there exists k ≥ 0 such that F(f)^k(S) = S′.

For any trap space T and its image T′ = F(f)(T), we have T′ ⊆ T, since by definition F preserves the fixed components of T. The free components in T might get fixed or remain free depending on the Boolean function f. In fact, T′ is a trap space of f, since for any fixed variable i ∈ V in T′, fi(x)=Ti′ by definition of F. Moreover, F(f)^k(T) is a trap space for any k ≥ 0. An example of a subspace Ω = 1∗∗∗ percolating to another subspace S = 10∗∗ is shown in Figure 1. Both Ω and S are trap spaces of f_↾Ω.

Remark 3.4. Let f be a Boolean function and S ⊆ 𝔹ⁿ a subspace. Let k ≥ 0 and Sk=F(f↾S)k(S). Then S^k is a trap space of f_↾S for every k ≥ 0.

Note that the paths in the dynamics of F(f) starting at a trap space T cannot have cycles and, consequently, all the reachable attractors from T in these dynamics are fixed points. When considering the synchronous dynamics of F(f), each initial trap space T leads to a unique fixed point.

Definition 3.5. Given a Boolean function f and a trap space T, we call the unique fixed point T′ of the synchronous dynamics of F(f) reachable from T the percolated subspace from T with respect to f, that is, T′ = F(f)^k(T) with k such that F(f)^k(T) = F(f)^r(T) for all r ≥ k.

In order to relate value percolation to control strategies, we derive some dynamical properties of percolated subspaces.

Lemma 3.6. Let f be a Boolean function, T ⊆ 𝔹ⁿ a trap space. Let k ≥ 0 and T^k = F(f)^k(T). Then for every x ∈ T there exists a path in D(f) from x to some y ∈ T^k.

Proof. It is enough to show that if T is a trap space, then for every x ∈ T there exists a path in D(f) from x to some y ∈ F(f)(T). Set T′ = F(f)(T), with I′ ⊆ V being the set of fixed variables in T′. Since T′ ⊆ T, for all i ∈ I′, fi(x)=Ti′ by definition of F. Now let us look at every update separately. In the asynchronous dynamics, for every i ∈ I′, x admits a successor y in AD(f) with yi=Ti′ and y_j = x_j for j ≠ i; therefore there exists a path from any state in T to T′. In the synchronous dynamics, fi(x)=Ti′ for all i ∈ I′ and so x admits a successor y ∈ T′. The case of the generalized asynchronous dynamics follows from the other cases, since all the paths in SD(f) or AD(f) are also paths in GD(f).

Corollary 3.6.1. Let f be a Boolean function, S ⊆ 𝔹ⁿ a subspace. Let k ≥ 0 and Sk=F(f↾S)k(S). Then for every x ∈ S there exists a path in D(f_↾S) from x to some y ∈ S^k.

Lemma 3.7. Let f be a Boolean function and S ⊆ 𝔹ⁿ a subspace. Let k ≥ 0 and Sk=F(f↾S)k(S). Then A is an attractor of f_↾S if and only if A⊆Sk and A is an attractor of f↾Sk.

Proof. As noted in Remark 3.4, S^k is a trap space of f_↾S. Then f↾S(x)=f↾Sk(x) for all x ∈ S^k. Therefore, any attractor of f↾Sk is also an attractor of f_↾S and if A is an attractor of f_↾S and A⊆Sk, then A is also an attractor of f↾Sk. Let A be an attractor of f_↾S. Then A=Reachf↾S(A). By Corollary 3.6.1, for every x ∈ S there exists a path in f_↾S from x to some y ∈ S^k. Therefore, Reachf↾S(A)∩Sk≠∅ and so, Sk∩A≠∅. Since S^k is a trap space of f_↾S, A⊆Sk. Therefore, all the attractors of f_↾S are contained in S^k and, consequently, are also attractors of f↾Sk.

Corollary 3.7.1. Let f be a Boolean function and S ⊆ 𝔹ⁿ a subspace and P ⊆ 𝔹ⁿ a subset. Let k ≥ 0 and Sk=F(f↾S)k(S). S is a control strategy for P if and only if S^k is a control strategy for P.

Corollary 3.7.1 provides a way to identify control strategies or discard candidate subspaces by using value percolation. Moreover, checking whether the percolated subspace satisfies the conditions of Definition 3.1 instead of the original subspace allows us to reduce the dimension of the restricted network and, consequently, to simplify the verification problem. In the example shown in Figure 1, it would be enough to check the subspace S to know whether the subspace Ω is a control strategy.

Percolation-only methods select candidate subspaces and percolate them. If the resulting subspace is contained in the target subspace, the candidate subspace is identified as a control strategy [10]. This type of control strategies can be identified efficiently [12]. However, additional control strategies might exist, that do not directly percolate to the target subspace, as shown in [9], or even to non-trivial intermediate trap spaces. An example of this scenario can be seen in Figure 2, where the control strategy Ω does not percolate to the target nor to any non-trivial trap space. With our approach, value percolation can also be exploited as a first step toward control strategy identification, as a means to achieve dimensionality reduction, as described by Corollary 3.7.1.

To improve control detection, we propose to define sufficient conditions on the system restricted to a candidate subspace to identify this candidate as a control strategy. Moreover, conditions on minimal trap spaces can be defined in order to deduce properties of the system attractors, in particular, their belonging to a target subset.

Definition 3.8. A set of trap spaces T is complete in D(f) if and only if for every attractor A of D(f) there exists T∈T such that A⊆T. A Boolean function f is complete in the dynamics D(f) if its minimal trap spaces are complete in D(f).

Completeness of the minimal trap spaces has been used for attractor approximation and it can be detected using model checking as described in [16]. The following proposition presents sufficient conditions for a subspace to be a control strategy. Given a candidate subspace, if the set of minimal trap spaces of the restricted system is complete and contained in the target subset, then the candidate subspace is a control strategy for that subset.

Proposition 3.9. Let f be a Boolean function, P, Θ ⊆ 𝔹ⁿ subspaces and T the set of minimal trap spaces of f_↾Θ. If all the trap spaces of T are contained in P and T is complete in D(f_↾Θ), then Θ is a control strategy of P.

Proof. Let A be an attractor of D(f_↾Θ). Since T is complete in D(f_↾Θ), there exists a minimal trap space T∈T such that A⊆T. Therefore, A⊆T⊆P.

Proposition 3.9 provides sufficient conditions that allow us to identify new control strategies missed by percolation-based approaches. An example of such a control strategy is shown in Figure 2. The subspace Ω is not a trap space nor percolates to any smaller subspace. Since f_↾Ω is complete in AD(f_↾Ω) and its only attractor is included in the target P, Ω is identified as a control strategy. We refer to this approach for control strategy identification as the completeness approach. However, it still does not characterize all the possible control strategies satisfying Definition 3.1. Figure 3 shows a subspace Ω that is a control strategy for the target P but, since Ω is the unique trap space in f_↾Ω and Ω ⊈ P, Proposition 3.9 cannot be applied. Thus, Ω is not detected as a control strategy by the completeness approach. To obtain the full solution set, we formulate a model checking approach, as shown in the next section.

This section provides a practical introduction to the model checking concepts required for the description of our approach. For a more extensive and detailed explanation of model checking we refer the reader to [17]. Model checking is a formal method used in computer science to solve verification problems. Its application to the control strategy problem presents many advantages, for instance the use of symbolic representation, which allows one to deal with systems with a large number of states, like STGs of Boolean networks. Moreover, many efficient algorithms have been developed and are available for running model checking queries. An overview of existing model checking tools in the context of biochemical networks analysis can be found in [15].

Model checking allows one to verify whether a given transition system satisfies a specific property. A transition system is defined as a set of states and a set of transitions, which represent changes from one state to another. Formally, a labeled transition system (LTS) is defined by a tuple (S,T,L) where S is a finite set of states, T ⊆ S × S is a transition relation such that (x¹, x²) ∈ T if there exists a possible transition from state x¹ to state x² and L : S → 2^AP is a labeling function with AP a finite set of atomic propositions. In the following, a transition (x¹, x²) will also be denoted by x¹ → x². The labeling function L gives a set L(x) ∈ 2^AP of atomic propositions for each state x which includes exactly the atomic propositions satisfied by x. In the Boolean context, an STG defines an LTS, where the set of states is 𝔹ⁿ and the transitions are defined by the Boolean function and the type of update that is chosen. For our purposes, we need a deadlock-free transition system, so we add extra transitions (x, x) ∈ T for every steady state x ∈ 𝔹ⁿ. We use the atomic propositions AP = {(v = c) | v ∈ V, c ∈ 𝔹} and define the labeling function by (v = c) ∈ L(x) if and only if x_v = c.

There are different ways to express properties of a transition system. In our case, we use Computational Tree Logic (CTL). CTL is based on a branching notion of time, where the behavior of the system is represented by a tree of states. In the case of Boolean networks, one can imagine that every path starting in a state x ∈ 𝔹ⁿ is represented as a branch in a tree rooted in x. In the following we introduce the main concepts of CTL.

We distinguish between state properties and path properties. In this context, a path is an infinite sequence x⁰, x¹, … ∈ S such that (xⁱ⁻¹, xⁱ) ∈ T for all i ≥ 1. A statement about a state or a path can be made using a CTL formula. A CTL state formula ϕ over the set of atomic propositions AP is of the form ϕ:=a |¬ϕ |ϕ1∨ϕ2 |ϕ1∧ϕ2 |Eφ |Aφ where a ∈ AP is an atomic proposition, E is the exists operator, A is the for all operator, ϕ, ϕ₁, and ϕ₂ are CTL state formulas and φ is a CTL path formula, which in our work will be of the form:

where F is the future operator, G the global operator, and ψ a CTL state formula. Note that we do not use the full expressiveness of CTL but only a subset of operators necessary to formulate our control queries. In accordance to the definition of semantics of CTL formulas, a state formula is evaluated for a state whereas a path formula is evaluated for an infinite path. When ϕ = a, where a is an atomic proposition, if a ∈ L(x) we say that the CTL state formula ϕ is satisfied at the state x and we write x ⊧ ϕ. For example, if ϕ = (i = 1) for some i ∈ V, then x ⊧ ϕ if (i = 1) ∈ L(x), that is x_i = 1. Analogously, we write π ⊧ φ when the CTL path formula ϕ is satisfied by a path π. See Table 1 for the recursive definition of the satisfaction relation ⊧ for transition systems and CTL formulas used in this work. Since we only use atomic propositions of the form (v = c) where v ∈ V and c ∈ 𝔹 and (v = c) ∈ L(x) if and only if x_v = c, atomic propositions can be interpreted in each state x according to the values of each variable in x. Thus, to simplify the notation, given an atomic proposition ϕ = (v = c), we define ϕ(x) = (x_v = c) so that ϕ is satisfied by x if and only if ϕ(x) = true. Figure 4 shows some examples of state and path formulas which are satisfied in an STG.

In this section, we present the basis of our new approach for the identification of all the minimal control strategies, based on model checking. To do so, we express the definition of control strategy in terms of CTL formulas. We start by rewriting it in terms of paths.

Lemma 4.1. Let f be a Boolean function, Θ ⊆ 𝔹ⁿ a subspace and P ⊆ 𝔹ⁿ a subset. The following are equivalent:

Θ is a control strategy for P in D(f).

Proof. (⇒) Let x ∈ Θ and let A be an attractor of D(f_↾Θ) that can be reached from x. Since Θ is a control strategy, A⊆P. Take y∈A. Since A is reached from x, there exists a path from x to y∈A⊆P and there are no paths from y leaving P.

(⇐) Let A be an attractor of f_↾Θ. Let x∈A. Since x∈A⊆Θ, there exists y ∈ P such that there exists a path in D(f_↾Θ) from x to y and there does not exist any path in D(f_↾Θ) from y to any state outside P. Since A is an attractor, y∈A and ReachD(f↾Θ)(y)=A. Then, A∩𝔹n\P=∅, that is, A⊆P, and Θ is a control strategy for P.

Before expressing Lemma 4.1 in terms of CTL formulas, we introduce a state formula ψ_S that is satisfied at a state x if and only if the state x belongs to the subspace S = Σ(I, c): ψS(x)=∧i∈I(xi=ci). This formulation can be extended to subsets as well. Clearly, every subset can be written as the union of subspaces, since a singleton set constitutes a subspace. Let P ⊆ 𝔹ⁿ be a subset. We define ϕ_P to be satisfied at a state x if and only if the state x belongs to x ∈ P:

ϕP(x)=∨S∈SψS(x) where S is a subspace cover of P.

Now we can express Lemma 4.1 in terms of CTL formulas, using ϕ_P(x) as defined above.

Lemma 4.2. Let f be a Boolean function, Θ ⊆ 𝔹ⁿ a subspace and P ⊆ 𝔹ⁿ a subset. The following are equivalent:

Θ is a control strategy for P in D(f).

Proof. Φ_P is satisfied at a state x if and only if there exists a path x = x⁰, x¹, … such that AG(ϕ_P) is satisfied at xⁱ for some i ≥ 1. Let y = xⁱ. AG(ϕ_P) is satisfied at y if and only if for all paths y = y⁰, y¹, …, for all i ≥ 0, ϕ_P is satisfied at yⁱ, that is, yⁱ ∈ P. Thus, by Lemma 4.1, Φ_P is satisfied for all x ∈ Θ if and only if Θ is a control strategy for P in D(f).

Note that Φ_P is not affected by the presence of non-attractive cycles since, by definition, they contain at least one state with an outgoing trajectory leading to an attractor.

The CTL formula Φ_P defined in Lemma 4.2 provides a way to determine whether a candidate subspace is a control strategy for P. The next section presents the implementation of this idea for control strategy identification.

Building on the model checking formulas derived in the previous section, we develop a method for control strategy identification. The formula derived in Lemma 4.2 can be used to define a CTL query that can determine whether a candidate subspace is a control strategy for a target subset. In addition, in order to improve the performance of the method, preliminary checks on the candidate subspace and the restricted network can be conducted to possibly discard it without exhaustive exploration. Moreover, the dimension of the problem, that is, the free variables of the Boolean function, can be reduced by restricting the function to the percolated subspace instead of the candidate subspace, as stated in Section 3.1. The complete implementation of the control strategy identification is detailed in Algorithm 1 and explained in the following.

The main algorithm takes as inputs a constant-free Boolean function f (i.e. all coordinate functions of f are non-constant), a target subset P and the type of update (asynchronous, synchronous or generalized) and returns the minimal control strategies for P in the respective dynamics. If the original Boolean function is not constant-free, a pre-processing step is applied in which the constant function values are percolated. Instead of the original function we then consider its restriction to the corresponding percolated subspace. As discussed in Section 3.1, this pre-processing does not impact the attractors of the system.

Since the aim is to identify optimal control strategies, the candidate subspaces S are taken randomly fixing an increasing number of variables, so that supersets of sets already defining a successful intervention are not considered (lines 5–8). Furthermore, an upper bound m for the size of the control strategies can be set. Moreover, the decisions made for each percolated subspace are stored in two variables ST (for positively checked subspaces) and SF (for negatively checked subspaces) to avoid repeating the same verification query.

The algorithm presented above is implemented using PyBoolNet [18], a Python package that allows generation and analysis of Boolean networks. PyBoolNet uses NuSMV to decide model checking queries for Boolean networks. It also provides an efficient computation of trap spaces for relatively large networks.

The network considered in this case study was recently introduced by Selvaggio et al. [20] to model how microenvironmental signals influence cancer-related phenotypes along the epithelial-to-mesenchymal transition (EMT). The original network consists of 56 components, ten of them being inputs and two readouts or outputs (see Figure 5). Since the original model is multivalued, we work with its booleanized version obtained with GINsim [19]. This booleanization maps a multivalued component of maximum value m to m Boolean components. For instance, a component taking the values 0, 1, 2, 3 is encoded using 3 Boolean variables that would take values 000, 100, 110, 111, respectively (see Table 2). Although this method introduces states that do not correspond to any value of the multivalued variable (non-admissible states), these cannot be part of any attractor since they always have at least one path to an admissible state and do not have incoming transitions from admissible states. Therefore, the asymptotic behavior generated strictly replicates the original model. The booleanized network of this case study consists of 60 Boolean variables, whose regulatory functions can be found in the PyBoolNet repository [18].

The asynchronous dynamics has 1,452 attractors, all of them steady states. They are classified according to the values of the readout components AJ and FA, which represent the different degrees of cell adhesions by adherens junctions and focal adhesions, respectively [20]. The eight resulting biological phenotypes are divided in four groups: epithelial (E1), mesenchymal (M1, M2, M3), hybrid (H1, H2, H3), and unknown (UN) (see Table 2).

To show the flexibility of our method, we analyse different control targets. We start by targeting single steady states (attractor control). Then, we target the subspaces corresponding to each phenotype (target control). Finally, we study the avoidance of the hybrid phenotype, setting as target the complement of the general hybrid phenotype.

In this section, we compare the model checking approach to other control methods currently available. We show that our method is able to capture all the minimal control strategies identified by other methods, uncovering in some cases minimal control strategies that might be missed by them.

In order to be able to compare different approaches, certain common features need to be chosen. Here, we consider control for any possible initial state. We separately compare to methods tackling attractor control and target control. Although approaches for target control can also be used for attractor control when the target attractor is a steady state or minimal trap space, they are usually aimed at targeting larger subspaces, determined for example by a phenotype, which often include several attractors. For this reason, we consider two different scenarios: one for attractor control and one for target control. The case of an arbitrary subset as target could not be considered for comparison, since no other method, to our knowledge, allows this possibility.

The comparison presented here encompasses each of the main approaches for control strategy identification discussed in previous sections (an overview of the main features of the control methods is shown in Table 4):

Since the methods for attractor control considered here (BA and SM) only work for asynchronous update, the comparison is only made for this dynamics. In the case of target control, control strategies are identified for both synchronous and asynchronous dynamics.

In view of the different nature of each method, we do not compare their running times. Some approaches require the computation of the system attractors to identify the control strategies (BA, SM). Others do not allow the choice of a certain attractor as target and identify control strategies for all the attractors simultaneously (SM). Therefore, a fair comparison with respect to the running times is hard to achieve. For this reason, we focus on the amount and size of the control strategies identified, provided that the program terminates within several hours. In particular, our method takes between a few minutes and 2 h for the networks and targets chosen for attractor control and between half an hour and 12 h for the networks and targets chosen for target control.

In order to capture different control scenarios, several biological networks of different sizes with different type and number of attractors are considered. A short description of each network is provided below. See Table 5 for an overview of the networks and their features. The Boolean rules for each biological network can be found in the PyBoolNet repository [18].

As mentioned in Section 4.3, our method takes as input a constant-free Boolean function. Therefore, all the networks with constant coordinate functions considered for comparison have been pre-processed so that the constant values are percolated and removed from the network. In such cases, for the sake of the comparison, all the methods have taken as input the reduced network.