Control Strategy Identification via Trap Spaces in Boolean Networks
CControl Strategy Identification via Trap Spacesin Boolean Networks
Laura Cifuentes Fontanals , Elisa Tonello , and Heike Siebert Freie Universit¨at Berlin, Germany Max Planck Institute for Molecular Genetics, Berlin, Germany
Abstract
The control of biological systems presents interesting applications such as cell reprogrammingor drug target identification. A common type of control strategy consists in a set of interventionsthat, by fixing the values of some variables, force the system to evolve to a desired state. This workpresents a new approach for finding control strategies in biological systems modeled by Booleannetworks. In this context, we explore the properties of trap spaces, subspaces of the state spacewhich the dynamics cannot leave. Trap spaces for biological networks can often be efficientlycomputed, and provide useful approximations of attraction basins. Our approach provides controlstrategies for a target phenotype that are based on interventions that allow the control to beeventually released. Moreover, our method can incorporate information about the attractors tofind new control strategies that would escape usual percolation-based methods. We show theapplicability of our approach to two cell fate decision models.
The control of biological systems presents interesting applications such as cell fate reprogramming,drug target identification for disease treatments or stem cells programming [5, 16]. Controlling acell fate decision network could for instance allow, in the case of cancer cells, to lead the system toan apoptotic state and, therefore, evolve towards the elimination of pathological cells [1]. Findingadequate candidates for control is a complex problem, in particular since the experimental testing ofall the possibilities is not feasible. Mathematical modeling can help address this problem by enabling in silico identification of possible effective candidates.Modeling of biological processes is often challenged by the lack of information about kinetic pa-rameters or specific reaction mechanisms. The Boolean formalism aims at capturing the qualitativebehavior of systems via a coarse representation of the relationship between the species of interest.Mechanisms underlying activation and inhibition processes are summarized in logical functions, allow-ing for two activity levels for each variable. The two values can represent for example if a gene isexpressed or not, or if the concentration of a protein is above or below a certain threshold. Booleanmodeling has in many instances been shown to capture the fundamental behaviors and dynamics ofbiological systems and has been widely used to make predictions or design strategies for therapeuticinterventions [3, 6, 7].Control of biological systems is a broad field that encompasses a variety of approaches and goals.Attractor control aims at leading the system to a desired attractor, starting from a particular initialstate (source-target control) [13] or from all possible initial states (full-network control) [19]. However,it is often useful to induce a desired phenotype rather than a specific attractor. Phenotypes are usuallydefined in terms of some biomarkers i.e., observable and measurable components that represent the1 a r X i v : . [ q - b i o . M N ] M a y ain characteristics of biological processes. The approach that focuses control on a set of relevantvariables is also known as target control [15, 17]. In this work, we are interested in full-network controlfor a target phenotype.There are different approaches for system interventions, that is, the way the control is applied tobiological systems. In the context of Boolean modeling, we consider as interventions the perturbationsor modifications that fix the value of some components (node control) [13, 19]. In the example of agene regulatory network, fixing a variable to a certain value can be understood as the knockout orpermanent activation of a gene. Among other approaches to Boolean network control is edge control,which targets the interactions between variables [2, 14]. For a gene regulatory network, edge controlcan be interpreted for instance as the modification of a protein to alter its interaction with a certaingene.Control of dynamical systems has been a popular research field in systems biology in the last years,also in the Boolean setting. Many approaches focus on the structure and topology of the network,for example by looking at feedback loops [18] or stable motifs [19], and several studies discuss thecomplexity and characteristics of such problems [8, 12]. Other approaches include techniques based ontopological information to reduce the size of the search space [15] or computational algebra methods[14]. Recent works have explored attractor control through the characterization of basins of attraction,that is, sets of states from which only a certain attractor can be reached [13]. However, the identificationof basins of attraction might require the exploration of the complete state space. Attractor reachabilitycan be investigated using trap spaces, which are subspaces that trajectories cannot leave. By definition,every trap space contains at least one attractor and, therefore, in some cases minimal trap spaces canbe good approximations for the attractors [10]. The identification of trap spaces in biological systemscan often be performed efficiently by exploiting properties of the prime implicants [9].Our approach aims to identify strategies for phenotype control by exploiting properties of trapspaces. We introduce the concept of space of attraction, a subspace that approximates the basinof attraction, to find control strategies without the need of computing the whole basin. We extendthis idea to define spaces of attraction for trap spaces and relate them to control strategies, which aredefined as sets of constraints that fix the value of some variables and induce a certain target phenotype.We exploit properties of trap spaces and computation techniques for target control to define a newmethod to compute control strategies that do not require a permanent intervention and allow thecontrol to be eventually released. Our approach can incorporate information about the attractorsto obtain new control strategies that might escape percolation-based target control techniques. Themethod presented here is widely applicable to Boolean models of biological systems and can provide,under certain conditions, control strategies that are independent of the type of update used in themodel.We start by giving a general overview about Boolean modeling (Section 2). Then we introduce theconcepts of control strategy and space of attraction in this setting (Section 3), providing the theoreticalbases for the computation of some types of control strategies. In Section 4, we present a method tocompute control strategies based on the theoretical principles explained in Section 3 and implementedusing the prime implicants of the function. Lastly, in Section 5 we show the applicability of our methodto two cell fate decision networks [7, 20]. A Boolean network on n variables is defined as a function f : B n → B n , where B = { , } . V = { , ..., n } is the set of variables of f , B n is the state space of the Boolean network and every x ∈ B n is a state ofthe state space. For any x ∈ B n and I ⊆ V , ¯ x I is defined as ¯ x Ii = x i for i ∈ V \ I and ¯ x Ii = 1 − x i for i ∈ I . If I = { i } , ¯ x I is written as ¯ x i .A dynamics on B n or state transition graph is a directed graph with vertex set B n . There are severalways of associating a dynamics to a Boolean network f . In the general asynchronous dynamics or2 eneral asynchronous state transition graph GD ( f ) there exists an edge from a vertex x to a vertex y ifand only if there exists ∅ (cid:54) = I ⊆ V such that ¯ x I = y and f i ( x ) = y i for every i ∈ I . Note that the generalasynchronous dynamics considers transitions which update subsets of components simultaneously ina non-deterministic way. By choosing different types of updates, other state transition graphs canbe defined. The asynchronous dynamics AD ( f ) is defined by considering the transitions updatingonly one component at a time and the synchronous dynamics SD ( f ) considers only the transitionswhere all the components that can be updated are updated at once. Note that AD ( f ) and SD ( f )are subgraphs of GD ( f ). To simplify the notation, D ( f ) will denote any of these dynamics associatedto f . The choice of asynchronous and general asynchronous updates is motivated by the attempt tocapture different, and sometimes unknown, time scales that might coexist in the modeled system. Anexample of asynchronous dynamics of a Boolean network is shown in Figure 1.A trap set T ⊆ B n is a set such that for all x ∈ T , if y is a successor of x in the dynamics, then y ∈ T . A minimal trap set under inclusion is an attractor . An attractor can be a stable state (or fixed point ), when it consists only of one state, or a cyclic (or complex) attractor when it is larger.In biological systems, stable states can be identified with different cell fates or cell types, and cyclicattractors with cell cycles or specific cell processes. Given a Boolean function f and an attractor A ,the weak basin of attraction of A is defined as the set of states x such that there exists a path from x to an element of A in D ( f ). The strong basin of attraction of A is the set of states in the weak basin of A that do not belong to the weak basin of attraction of any other attractor different from A . Figure 1shows the weak and strong basins for an attractor in an asynchronous state transition graph.The control interventions considered in this work consist in fixing the values of some components.Formally, given a state c ∈ B n and a subset of variables I ⊆ V , we define the subspace induced by c and I as the set Σ( I, c ) = { x ∈ B n | ∀ i ∈ I, x i = c i } . The variables in I are called fixed variables ,while the other variables are called free . We denote subspaces as states, using the symbol ∗ for thefree variables. For example, the subspace { x ∈ B | x = 1 and x = 0 } is denoted as 1 ∗ ∗ .The identification of control variables requires examining the effect that fixing certain variableshas on the dynamics. Given a Boolean function f and a subspace Θ = Σ( I, c ), the restriction of thefunction f to the subspace Θ is defined as: f (cid:22) Θ : Θ → Θ , where for all i ∈ V , ( f (cid:22) Θ ) i ( x ) = (cid:26) f i ( x ) , i / ∈ I,c i , i ∈ I. Note that f (cid:22) Θ : Θ → Θ can be identified with a Boolean network g : B m → B m , where m = n − | I | . Viathis identification, we extend all the definitions that apply to a Boolean network to such restrictions.For example, the state transition graph corresponding to f (cid:22) Θ : Θ → Θ is defined as usual, only withvertex set Θ instead of B n (see Figure 2). Moreover, if T is a trap set in D ( f ), then T ∩ Θ is a trapset in D ( f (cid:22) Θ ).A subspace that is also a trap set is called a trap space . While trap sets and attractors might varywhen considering different types of dynamics, trap spaces are independent of the type of update. TheBoolean function represented in Figure 1 has four trap spaces: 000, 111, 0 ∗ ∗ ∗ ∗ .In this work we aim at using trap spaces to find control strategies for phenotypes. Phenotypesare usually defined in terms of the state of some measurable components called biomarkers, which areobservable components that can be used as indicators of different cell types or cell fates or to distinguishbetween healthy and pathological conditions. Although the notion of phenotype is usually related tostability, we extend this concept to consider any possible state in order to allow non-attractive statessatisfying the phenotype characteristics to become attractors in the controlled system. Thus, in thiswork, we define a phenotype as a subspace. 310 111100 101010 011000 001 Basins of attraction of A : • Strong ( A ) = { , , , , }• W eak ( A ) = { , , , , , , } Spaces of attraction of A : • Ω = 0 ∗ ∗ , Ω = 00 ∗ , Ω = 01 ∗ , Ω = 0 ∗
0, Ω = 0 ∗
1, Ω = ∗ = 000, Ω = 001 , Ω = 010, Ω = 011, Ω = 101, withΩ i (cid:40) Strong ( A ) for all 1 ≤ i ≤ Figure 1:
Asynchronous dynamics of the Boolean function f ( x ) = (¯ x ¯ x x ∨ x x , x ¯ x ¯ x ∨ x x x , x x ∨ x x ∨ x x ), with attractors A = 000 and A = 111 and trap spaces 000, 111, 0 ∗ ∗ ∗ ∗ . Allthe spaces of attraction of A are included in its strong basin (in red) while the basin itself is not aspace of attraction.
000 001 100 110000 010
Figure 2:
Asynchronous dynamics of the Boolean function f ( x ) = ( x ¯ x ∨ ¯ x ¯ x , x ∨ x , x x ∨ x x )(left) and f (cid:22) Ω ( x ) = ( x ∨ ¯ x , x , 0) with Ω = ∗ ∗ P = { } in AD ( f ). Ω does not percolate to P . The strong basin of attraction of an attractor A can be naturally related to control since, by definition,it contains all the states that have paths to A but not to any other attractor. In contrast to methodsrequiring basin exploration, we use subspace approximation of the basins combined with trap spacescomputation. To do so, we extend the notion of basin of attraction to trap sets. We then exploituseful properties of trap spaces, e.g. independence of the update, efficient identification and potentialapproximation of attractors, to develop a new approach for the identification of control strategies. We now formalise the notion of control strategy for a phenotype. A control strategy is a subspacedefined by a set of interventions that fix the value of some variables and thus force all attractors to becontained in the subspace defining the phenotype.
Definition 3.1.
Given a Boolean function f and a phenotype P ⊆ B n , a control strategy (CS) for thephenotype P in D ( f ) is a subspace Θ ⊆ B n such that, for any attractor A of D ( f (cid:22) Θ ), A ⊆ P .If the desired phenotype is a stable state in the original dynamics ( P = { y } , y ∈ B n ), a controlstrategy for P is a subspace Θ such that y is the only attractor of f (cid:22) Θ . Figure 2 shows an example ofa control strategy for a stable state. The size of the subspace defining a control strategy representsthe number of interventions in the system. Therefore, the most interesting control strategies are thesubspaces that are maximal with respect to inclusion.4 common approach in the context of control is the use of value percolation [15, 17]. Differentcombinations of variables to be fixed are considered, and their values propagated iteratively until aninvariant subspace is reached. A combination of variables and values is an intervention strategy if thesubspace obtained at the end of the iterative percolation process is contained in the target phenotype.Strategies obtained with this approach satisfy the conditions of Definition 3.1. However, the class ofcontrol strategies identified by the definition is larger, as we will discuss in the following. Trap sets are sets of states that the dynamics cannot leave. Each trap set contains, as a consequence,at least one attractor. The concept of basin of attraction defined for an attractor can be naturallyextended to trap sets. As mentioned before, we wish to approximate basins of attraction by subspaces.Combining these two ideas, we introduce the concept of space of attraction of a trap set T as a subspaceΩ such that from any state in Ω there exists a path to T and no trap set disjoint from T is reachablefrom Ω. Definition 3.2.
Let f be a Boolean function and T be a trap set of f . A space of attraction of thetrap set T in D ( f ) is a subspace Ω such that for all x ∈ Ω and for any trap set S , if there exists a pathin D ( f ) from x to an element of S , then S ∩ T (cid:54) = ∅ .Definition 3.2 implies the existence of a path from the space of attraction Ω to the trap set T . Atrap set can have many spaces of attraction. In fact, any subspace contained in a space of attractionis also a space of attraction. Moreover, if there is only a unique trap set T m minimal with respect toinclusion contained in a trap set T , any space of attraction of T is also a space of attraction of T m .Both trap spaces and spaces of attraction are subspaces that characterize the long term behavior ofthe system. However, in contrast to trap spaces, spaces of attraction can depend on the update.If a trap set is an attractor, its spaces of attraction can be related to its basins of attraction.The spaces of attraction of an attractor A are clearly contained in the strong basin of A since, byDefinition 3.2, none of the other attractors can be reached from any state inside the space of attraction.However, the strong basin of attraction of A might not be a space of attraction (see Figure 1).Spaces of attraction, as well as basins, might include paths crossing non-attractive cycles in thestate transition graph. As a consequence, some paths starting in the space of attraction (or basin)might not reach the trap space (or attractor), staying indefinitely in non-attractive cycles. While invery specific circumstances such behavior might be relevant, generally it constitutes an artifact arisingfrom the non-deterministic update. Here, we extend the standard view on basins of attraction tospaces of attraction, assuming the trajectories of interest will eventually leave non-attractive stronglyconnected components in the state transition graph.The condition that a subspace needs to satisfy to be a space of attraction of a trap set T getssimplified when the subspace considered is the entire state space. In this case, it is only required thatfrom every state in the state space trap set T can be reached (see Lemma 3.3), since it immediatelyfollows that there cannot be a trap set disjoint from T . Lemma 3.3.
Let f be a Boolean function and T a trap set of f . Then B n is a space of attraction ofthe trap set T in D ( f ) if and only if for all x ∈ B n there exists a path in D ( f ) from x to some y ∈ T . The application of Lemma 3.3 to the restriction on a subspace immediately yields the followingcorollary.
Corollary 3.3.1.
Let f be a Boolean function, T a trap set of f and Ω a subspace such that T ⊆ Ω .Then Ω is a space of attraction of T in D ( f (cid:22) Ω ) if and only if for all x ∈ Ω there exists a path in D ( f (cid:22) Ω ) from x to some y ∈ T . a)
000 001 100
000 010 (b)
000 001 100
000 010
Figure 3: (a)
Ω = ∗∗ AD ( f ) and AD ( f (cid:22) Ω ), with f ( x ) = (¯ x ∨ x ¯ x , x ¯ x ∨ x x , ¯ x ¯ x ∨ x x ) and f (cid:22) Ω ( x ) = (¯ x ∨ x , x , (b) Ω = ∗ ∗ AD ( g ) butnot for AD ( g (cid:22) Ω ), with g ( x ) = ( x ¯ x ∨ ¯ x x ∨ x ¯ x , x ¯ x ∨ x x , ¯ x ¯ x ∨ x x ) and g (cid:22) Ω ( x ) = ( x , x , f to lead thedynamics to a certain trap set. If T is a trap space, there is always a trivial space of attraction for therestricted function which is T itself.Note that a subspace Ω that is a space of attraction of T for the Boolean function f is not necessarilya space of attraction for the restricted function f (cid:22) Ω (see Figure 3).Given a trap space T that only contains attractors belonging to a certain phenotype P , any spaceof attraction that leads the system to T would also lead it to an attractor belonging to P . In otherwords, any space of attraction for a trap space T is also a control strategy for a phenotype P if T onlycontains attractors belonging to P . The following proposition formalizes this idea. Proposition 3.4.
Let P ⊆ B n be a subspace and f a Boolean function. Let T be a trap space suchthat if A ⊆ T is an attractor of D ( f ) , then A ⊆ P . Let Ω be a space of attraction of T in D ( f (cid:22) Ω ) suchthat T ⊆ Ω . Then Ω defines a control strategy in D ( f ) for P .Proof. Let A be an attractor for D ( f (cid:22) Ω ). Then A ⊆ Ω. Since Ω is a space of attraction of T in D ( f (cid:22) Ω )and A is a trap set in D ( f (cid:22) Ω ), T ∩ A (cid:54) = ∅ . As T and A are trap sets, T ∩ A is also a trap set in D ( f (cid:22) Ω ).Since A is minimal, A = T ∩ A ⊆ T . Then, since T is a trap space and for all x ∈ T, f (cid:22) Ω ( x ) = f ( x ), A is also an attractor of D ( f ) and, therefore, A ⊆ P .Since a trap space is always a space of attraction of itself, given a subset P ⊆ B n , any trap space T containing only attractors in P is a control strategy for P .The type of control strategies identified by Proposition 3.4 allow the interventions to be releasedafter a certain number of steps. That is because these control strategies induce the target phenotypeby leading the system to a trap space. Once the trap space is reached, since the dynamics cannot leaveit, the control can be released and the system will remain in the trap space, eventually evolving to thephenotype of interest. As explained in the previous section, control strategies for a phenotype P can be found by identifyingspaces of attraction of trap spaces containing only attractors in P . In this section, we explore ways offinding spaces of attraction for trap spaces.Given a trap space T , we look for a subspace Ω such that from all states in Ω there is a path to T in D ( f (cid:22) Ω ). To do so, we use the idea of value percolation, which is a common approach in the contextof control. As explained in Section 3.1, it is based on the fact that the constraints given by the fixedvariables of a subspace might induce further variables to get fixed. Thus, in our setting, a subspaceΩ = Σ( W, c ) that percolates to the trap space T = Σ( U, c ) is a space of attraction of T in f (cid:22) Ω . Thefollowing lemma formalizes this idea. 6 emma 3.5. Let f : B n → B n be a Boolean function, c ∈ B n and S = Σ( U, c ) , Ω = Σ(
W, c ) subspacesof B n such that S ⊆ Ω and W ⊆ U ⊆ V . If for all s ∈ U \ W , f s ( x ) = c s for all x ∈ Ω , then for all x ∈ Ω there exists a path in D ( f (cid:22) Ω ) from x to some y ∈ S .Proof. Since the proof depends on the update, we treat each case separately. D = AD : For each x ∈ Ω and for each s ∈ U \ W such that x s (cid:54) = c s , f s ( x ) = c s . Therefore, x admits a successor y in AD ( f (cid:22) Ω ) with y s = c s . This implies the existence of a path in AD ( f (cid:22) Ω ) fromany state in Ω to S . D = SD : For each x ∈ Ω and for each s ∈ U , f s ( x ) = c s . Therefore, x admits a successor y ∈ S ⊆ Ω in SD ( f (cid:22) Ω ). D = GD : Since all the paths in AD ( f ) and SD ( f ) are also paths in GD ( f ), the conclusion followsfrom the previous cases.Lemma 3.5 can be extended with Corollary 3.3.1 to provide conditions that allow the identificationof spaces of attraction. Lemma 3.6.
Let f : B n → B n be a Boolean function and T = Σ( U, c ) a trap space of f with U ⊆ V and c ∈ B n . Let Ω = Σ(
W, c ) be a subspace of B n such that T ⊆ Ω and W ⊆ U ⊆ V . If f s ( x ) = c s forall x ∈ Ω and s ∈ U \ W , then Ω is a space of attraction of T for D ( f (cid:22) Ω ) . To improve the spaces of attraction obtained with Proposition 3.4, we can extend Lemma 3.6applying the idea used in Lemma 3.5 several times, building a path of percolated subspaces ending inthe trap space T . Proposition 3.7.
Let f : B n → B n be a Boolean function and let c ∈ B n . Let T = Σ( U, c ) bea trap space and Ω = Σ(
W, c ) be a subspace containing T with W ⊆ U ⊆ V . Let I = W and I k +1 = { s ∈ U | s ∈ I k or f s ( x ) = c s for all x ∈ S k } , where S k = Σ( I k , c ) . If there exists a k T suchthat I k T = U , then Ω is a space of attraction of T for D ( f (cid:22) Ω ) . Proposition 3.7 gives sufficient conditions for a subspace to be a space of attraction of a trap spacein the restriction and, together with Proposition 3.4, provides a way to identify control strategies fora given phenotype. However, not all spaces of attraction fall under the conditions given by Proposi-tion 3.7. The example in Figure 3 (a) shows a space of attraction Ω = ∗ ∗ T = 110,which is also a control strategy for P = { } , where Ω does not percolate to T .Sometimes the attractors of a system of interest are known. In other cases they are not knownbut can be approximated by minimal trap spaces [10], that is, each minimal trap space containsonly one attractor and every attractor is included in a minimal trap space. This information is notusually exploited by target control methods, which often rely solely on percolation-like techniques.The approach described in this work can use this knowledge to find additional control strategies. Ifthe attractors are known or they can be approximated by minimal trap spaces, we can easily findtrap spaces satisfying the conditions of Proposition 3.4 by simply checking whether these attractorsor minimal trap spaces are included in a trap space. Therefore, larger trap spaces containing onlyattractors of the target phenotype can be identified. By Proposition 3.4, spaces of attraction for thesetrap spaces are also control strategies for the phenotype. These control strategies do not necessarilypercolate to the phenotype and, therefore, might not be identified by usual percolation techniques.Figure 2 shows an example of such a control strategy, where Ω = T = ∗ ∗ T , which contains only the attractor A = 110, and so, is a control strategy for thephenotype P = A . Note that Ω does not percolate to A .The attractors of a Boolean network might vary in different dynamics. Therefore, the trap spacessatisfying Proposition 3.4 and the control strategies characterized by them might also be dependenton the dynamics. Conversely, the spaces of attraction obtained by Proposition 3.7 are independentof the update. Thus, if the trap spaces considered satisfy the conditions of Proposition 3.4 in all thedynamics, the control strategies identified are also independent of the update.7 Computation of control strategies
We propose a method to find control strategies for a given phenotype, using the ideas explained in theprevious section. The main steps of the method are represented in Figure 4 and the detailed procedureis shown in Algorithm 1.In order to implement the computation of the control strategies, we use the prime implicants ofthe function. Given a Boolean function f : B n → B n , a c -implicant of f i , with c ∈ B and i ∈ V , is asubspace Q such that f i ( x ) = c for all x ∈ Q . A prime implicant is an implicant that is maximal underinclusion. Given T = Σ( U, c ), finding a subspace satisfying the hypothesis of Lemma 3.6 is equivalentto finding a subspace that is a c i -implicant of f i for all i ∈ U . Moreover, prime implicants can also beused to compute the trap spaces [9]. The computation of the prime implicants of a Boolean functionis in general a hard problem. However, networks modeling biological systems are usually relativelysparse, since the number of components regulating a variable is relatively small compared to the size ofthe network. Therefore, they are rather tractable in terms of prime implicants computation. Severaltools are available for the computation of prime implicants and trap spaces of Boolean functions. Weuse PyBoolNet [11], a Python package that allows generation and analysis of Boolean networks andprovides an efficient computation of prime implicants and trap spaces for quite large networks.
PhenotypeBoolean network Prime implicantsand trap spacesSelectedtrap spaces Spaces ofattractionControl strategies
Figure 4:
Main steps of the method for finding control strategies for a phenotype, represented incolor boxes according to their role: inputs (blue), precomputation (green), main computation (beige),output (red).We describe now the main steps of the method, outlined in Figure 4.
Algorithm 1
Control strategies for a phenotype P
Input : f Boolean function, P phenotype, attr attractors of f (optional) m limit size of the control strategies (optional) Output : control strategies for P function ControlStrategies ( f , P , attr ) T ← trapSpaces( f ) selTS ← selectedTrapSpaces1( T , P ) if attr (cid:54) = ∅ then : selTS ← selTS + selectedTrapSpaces2( T , P , attr ) CA ← ∅ for i in { , . . . , min( m , n ) } do : (cid:46) n total number of variables S ← { S subspace : | fixed(S) | = i , ∃ T ∈ selTS with T ⊆ S } for S in
S do : if (S (cid:54)⊆ S’ for all S’ in CA ) and isSpaceAttraction( f , S, selTS ) then : add S to CA return CA lgorithm 2 Subspace is a space of attraction
Input : f Boolean function, S subspace, TS trap spaces Output : True if S is space of attraction of a trap space in TS . False otherwise. function IsSpaceAttraction (f, S, TS ) f’ ← percolateFunction(f, S) return isNotEmpty( { T in TS : T ⊆ S and fixed(T) ⊆ fixed(f’) } ) Inputs.
The inputs are the Boolean function describing the system and the subspace of the targetphenotype P . The attractors, if known, are also used as input. Prime implicants and trap spaces canbe given as input or computed from the Boolean function. Selection of trap spaces.
Trap spaces of interest are divided into two types: trap spaces containedin P (Type 1) and trap spaces not contained in P but containing only attractors in P (Type 2). Astrap spaces have been identified in the previous step, this selection only requires checking whether atrap space belongs to one of the types (Algorithm 1: 3-5). Trap spaces of Type 2 are only identifiedwhen all the attractors are known or can be approximated by minimal trap spaces. In order to avoidunnecessary calculations, we do not consider trap spaces that percolate to smaller ones, since if a trapspace T percolates to a trap space T , all spaces of attractions of T are also spaces of attraction of T . Computation of spaces of attraction.
Spaces of attraction for the trap spaces from the previousstep are computed using the theoretical principles described in Proposition 3.7. The detailed procedureis shown in Algorithm 1: 6-11. For each subspace S that contains at least one of the selected trapspaces (Algorithm 1: 8), it is checked whether it is a space of attraction for one of the selected trapspaces (Algorithm 1: 10). To do so, the percolated function of f obtained by fixing the variables in S is calculated (Algorithm 2: 2). If T is contained in the subspace generated by S and all the variablesfixed in T are also fixed in the percolated function, then the subspace generated by S is a space ofattraction of T (Algorithm 2: 3). Since the aim is to find maximal spaces of attraction satisfyingthis property, the subspaces S are taken randomly fixing an increasing number of variables, so thatsupersets of sets already defining a space of attraction are not considered (Algorithm 1: 8, 10). Output
The obtained spaces of attraction are control strategies for the phenotype P by Proposi-tion 3.4 and, therefore, are returned as output.The method also allows to include some constraints on the control strategies. One example isthe exclusion of some components, which can be taken into account when selecting the subspaces S (Algorithm 1: 8). Another constraint to consider is on the size of the control strategies. Imposinga limit on the number of interventions might allow to reduce the computational cost without losinginteresting solutions, since small control strategies are usually the most relevant. In this section we discuss the application of our method to two Boolean networks describing cell fatedecision processes. In the first case study we consider two different control problems, one having aphenotype as target for the control, the second targeting single attractors. The second case studyfocuses on phenotype control. We compare the control strategies identified by our approach to theones obtained using exclusively value percolation, as described in Section 3.1. We show that, for bothexamples, new control strategies can be identified with the procedure introduced in this work.All computations in this section were done on an 8-processor computer, Intel(R)Core(TM) i7-2600CPU at 3.40GHz, 16GB memory, without any use of parallelization.9 .1 MAPK network
The network considered in this case study was introduced by Grieco et al. (2013) [7] to model theeffect of the Mitogen-Activated Protein Kinase (MAPK) pathway on cell fate decisions taken in patho-logical cells (see Figure 5). It uses 53 Boolean variables, four being inputs (DNA-damage, EGFR-stimulus, FGFR3-stimulus and TGFBR-stimulus) and three outputs (Apoptosis, Proliferation andGrowth-Arrest).The asynchronous dynamics has 18 attractors, 12 being stable states and 6 cyclic attractors. All ofthem can be approximated by minimal trap spaces, since each minimal trap space only contains oneattractor and there is no attractor that is not contained in a minimal trap space [10]. Therefore, wecan use trap spaces of both Type 1 and Type 2 to compute control strategies.The phenotype chosen as target for the control is the apoptosis phenotype, which is defined in [7]as the states fixing Apoptosis and Growth Arrest to 1 and Proliferation to 0. There are 103 non-percolating trap spaces, which are trap spaces that do not percolate to smaller ones, containing onlyattractors in the apoptosis phenotype. Of these, 64 are of Type 1 and 39 of Type 2. We set anupper bound of four components to the size of the control strategies, since generally only small controlstrategies are of interest and this limit already allows to find relevant ones. In addition, we excludeinterventions that fix any of the output nodes of the network. In this setting, we identify two controlstrategies of size 1 ( { TGFBR-stimulus = 1 } and { DNA-damage = 1 } ) and no control strategies of size2, 3 and 4. The running time is around 13 minutes.Using exclusively the percolation of the fixed values we identify two control strategies of size 1( { TGFBR-stimulus = 1 } and { TGFBR = 1 } ), 121 control strategies of size 2, 164 of size 3 and 139 ofsize 4. Looking at the Boolean function, we observe that TGFBR is uniquely regulated by TGFBR-stimulus, so fixing TGFBR-stimulus to 1 implies that TGFBR is also fixed to 1 and, therefore, theseinterventions are equivalent in terms of their effect on the apoptosis phenotype. However, it is obviousthat if the control fixing TGFBR to 1 is released, TGFBR could be updated to zero again by TGFBR-stimulus, and this change would induce the system to leave the apoptosis phenotype. Therefore, thecontrol of TGFBR requires a permanent intervention.Our method uncovers the control strategy { DNA-damage = 1 } , which is not obtained by usingsolely value percolation. In fact, the percolation of the subspace defined by this strategy does notreach the phenotype, but stops at the subspace T = { DNA-damage = 1, ATM = 1, TAOK = 1 } .However, since T is one of the trap spaces selected by our method, the constraint { DNA-damage = 1 } is identified as a control strategy.Of the control strategies of size 2, 3 and 4 that can be identified by percolation, 18, 13 and 7respectively are supersets of the control strategy { DNA-damage = 1 } identified by our method. Forthis reason, the subspaces obtained by percolating these interventions are contained in the trap space T mentioned above and therefore the associated control can be eventually released, without affectingthe reachability of the target. The remaining control strategies are not guaranteed to lead to a trapspace. As a consequence, in these cases, an early release of the control could lead to the loss of thecontrol goal. This illustrates how our method can complement previous approaches, by identifyingcontrol strategies of reduced complexity, and, consequently, reducing the number of interventions tobe considered, while at the same time providing information about the effects of a possible release ofthe control.The components appearing in the minimal control strategies identified (DNA-damage and TGFBR-stimulus) correspond to two inputs of the model. These inputs represent anti-proliferative stumuli fromthe MAPK network [7] and, therefore, can be expected to play an important role in the phenotypedecision. It is, however, certainly interesting that they are capable of fully inducing the apoptosisphenotype without further conditions on internal processes.In addition to the control problem for the apoptosis phenotype, we also searched for control strate-gies for the 10 apoptotic stable states. We set the maximum size of control strategies to five. Foreight stable states ( A to A in Table 1) exactly one control strategy of size 4 is obtained. For stable10tate A , two control strategies of size 5 are found, and for A no control strategies up to size 5 areidentified. The list of stable states and their control strategies can be found in the supplementarymaterial. The running time for one stable state is around 26 minutes.Since the chosen stable states belong to the apoptosis phenotype, all the selected trap spaces arealso considered when computing the control strategies for the apoptosis phenotype. Therefore, thecontrol strategies of the stable states are subspaces of the ones obtained for the apoptosis phenotype.One of the main differences is that the four inputs are present in all the control strategies of the stablestates. The input variables are, by definition, not regulated by any component, and therefore must bedirectly controlled if the value in a given steady state is to be achieved. The analysis of the controlproblem for the phenotype revealed that fixing DNA-damage to 1 is enough to lead the system tothe apoptosis subspace, but fixing the additional inputs is necessary to obtain a specific steady state.Fixing the four inputs is already enough to induce the stable states A to A solely by percolation.However, the stable states A and A require additional internal processes to be controlled. For A ,the two control strategies identified do not percolate directly to the attractor, but lead the dynamicsto one of the selected trap spaces. For A , no control strategies up to size 5 are found neither byour method nor percolation techniques, suggesting that a higher number of interventions might benecessary. These observations show that control for a phenotype can be more achievable than for aspecific attractor, and thus in some cases more interesting for application. We now consider a control problem for the network introduced by Zhang et al. (2008) [20] to modelthe T cell large granular lymphocite (T-LGL) survival signaling network (see Figure 6). It consistsof 60 Boolean variables, six being inputs (CD45, IL15, PDGF, Stimuli, Stimuli2 and TAX) and threereadouts (Apoptosis, Proliferation and Cytoskeleton-signaling).The asynchronous dynamics has 156 attractors, 86 being stable states and 70 cyclic attractors. Asin the previous network, all of them can be approximated by minimal trap spaces [10]. Therefore, wecan use trap spaces of both Type 1 and Type 2 to compute control strategies.We consider the apoptosis phenotype defined by fixing Apoptosis to 1 and Proliferation to 0. Notethat the third readout, Cytoskeleton signaling, is forced to 0 by its regulator Apoptosis having value1. There are 883 non-percolating trap spaces, which are trap spaces that do not percolate to smallerones, containing only attractors in the apoptosis phenotype. 729 trap spaces are of Type 1 and 154 ofType 2. As in the previous case study, we set an upper bound of four components to the size of thecontrol strategies and we exclude interventions that fix any of the readout nodes of the network. Inthis setting, six control strategies are identified: three of size 3 ( { CD45 = 0, IL15 = 0, PDGF = 1 } , { CD45 = 0, IL15 = 0, Stimuli = 1 } , { CD45 = 0, IL15 = 0, TAX = 1 } ) and three of size 4 ( { CD45 =1, PDGF = 0, PDGFR = 0, Stimuli2 = 1 } , { CD45 = 1, PDGF = 0, S1P = 0, Stimuli2 = 1 } , { CD45= 1, PDGF = 0, SPHK1 = 0, Stimuli2 = 1 } ). The running time is around 15 minutes.The three control strategies of size 3 consist only of input components. All the control strategiesof size 4 have three components in common while the fourth varies within PDGFR, S1P and SPHK1,suggesting that these three interventions might be equivalent in terms of their effect on the apoptosisphenotype. In fact, by looking at the Boolean function, we observe that fixing PDGFR = 0, impliesSPHK1 = 0, which also implies S1P = 0. Identifying such equivalent interventions a priori might allowto reduce the computational cost of the method.Using only percolation we find exactly one control strategy of size 1 ( { Caspase = 1 } ) and none ofsize 2, 3 or 4. However, this control strategy is relatively trivial since the Caspase component is directlyregulating Apoptosis. The control strategies identified by our method do not percolate directly to thephenotype. At the end of the percolation process, the dynamics reaches one of the trap spaces selectedas containing only attractors in the apoptosis phenotype.11his case study highlights the added value of our approach which can uncover relevant systeminterventions that are not identified by usual percolation approaches. In this work, we considered properties of trap spaces and principles of target control to introduce anew approach to compute control strategies. The procedure proposed is applicable to both phenotypeand attractor control and allows the interventions to be released after a certain amount of time, incontrast to usual target control methods that require permanent interventions.The approach presented here is widely applicable to Boolean models of biological systems and canprovide intervention strategies that are independent of the type of update considered in the modeling.Moreover, restrictions on the control strategies, in the form of variables to be excluded, can be added.Our approach also allows to incorporate information about the attractors, with the possibility toobtain control strategies that escape regular percolation-based techniques. As demonstrated with thetwo case studies, our method can identify new control strategies that require a small number of controlvariables, and can thus reveal valuable intervention approaches.Our approach efficiently identifies control strategies for relatively large biological networks. A nat-urally important further step is a rigorous comparison with existing methods, for instance approachesbased on stable motifs [17, 19]. Furthermore, the performance of the method could benefit from theadoption of fine-tuning strategies developed to speed up some of the procedures involved in candidatescreening. For instance, we could consider the reduction of the size of the search space by identifyinga priori equivalent interventions, adapting existing approaches [15]. Further steps also include theextension of the method to other types of control, such as edge interventions or sequential control.
Acknowledgments
E.T. was funded by the Volkswagen Stiftung (Volkswagen Foundation) under the funding initiativeLife? - A fresh scientific approach to the basic principles of life (project ID: 93063).
References [1] Baig, S., Seevasant, I., Mohamad, J., Mukheem, A., Huri, H.Z., Kamarul, T.: Potential of apop-totic pathway-targeted cancer therapeutic research: Where do we stand? Cell Death & Disease (1), e2850 (2016). https://doi.org/10.1038/cddis.2015.275[2] Biane, C., Delaplace, F.: Causal reasoning on boolean control networks based on abduction: The-ory and application to cancer drug discovery. IEEE/ACM Transactions on Computational Biologyand Bioinformatics (5), 1574–1585 (2019). https://doi.org/10.1109/TCBB.2018.2889102[3] Calzone, L., Tournier, L., Fourquet, S., Thieffry, D., Zhivotovsky, B., Barillot, E., Zinovyev, A.:Mathematical modelling of cell-fate decision in response to death receptor engagement. PLOSComputational Biology (3), 1–15 (2010). https://doi.org/10.1371/journal.pcbi.1000702[4] Chaouiya, C., Naldi, A., Thieffry, D.: Logical Modelling of Gene Regulatory Networks withGINsim., vol. 804, pp. 463–79 (2012) 125] Csermely, P., Korcsm´aros, T., Kiss, H.J., London, G., Nussinov, R.: Structure and dynamics ofmolecular networks: A novel paradigm of drug discovery: A comprehensive review. Pharmacology& Therapeutics (3), 333 – 408 (2013). https://doi.org/10.1016/j.pharmthera.2013.01.016[6] Flobak, ˚A., Baudot, A., Remy, E., Thommesen, L., Thieffry, D., Kuiper, M., Lægreid, A.: Discov-ery of drug synergies in gastric cancer cells predicted by logical modeling. PLOS ComputationalBiology (8), 1–20 (2015). https://doi.org/10.1371/journal.pcbi.1004426[7] Grieco, L., Calzone, L., Bernard-Pierrot, I., Radvanyi, F., Kahn-Perl`es, B., Thieffry, D.: Integra-tive modelling of the influence of mapk network on cancer cell fate decision. PLOS ComputationalBiology (10), 1–15 (10 2013). https://doi.org/10.1371/journal.pcbi.1003286[8] Kim, J., Park, S.M., Cho, K.H.: Discovery of a kernel for controlling biomolecular regulatorynetworks. Scientific Reports , 2223 (2013). https://doi.org/10.1038/srep02223[9] Klarner, H., Bockmayr, A., Siebert, H.: Computing maximal and minimal trap spaces of booleannetworks. Natural Computing , 535–544 (2015). https://doi.org/10.1007/s11047-015-9520-7[10] Klarner, H., Siebert, H.: Approximating attractors of boolean networks by itera-tive ctl model checking. Frontiers in Bioengineering and Biotechnology , 130 (2015).https://doi.org/10.3389/fbioe.2015.00130[11] Klarner, H., Streck, A., Siebert, H.: PyBoolNet: a python package for the genera-tion, analysis and visualization of boolean networks. Bioinformatics (5), 770–772 (2016).https://doi.org/10.1093/bioinformatics/btw682[12] Liu, Y.Y., Slotine, J.J., Barab´asi, A.L.: Controllability of complex networks. Nature , 167–173(2011). https://doi.org/10.1038/nature10011[13] Mandon, H., Su, C., Haar, S., Pang, J., Paulev´e, L.: Sequential reprogramming of booleannetworks made practical. In: Bortolussi, L., Sanguinetti, G. (eds.) Computational Methodsin Systems Biology. vol. 11773, pp. 3–19. Springer International Publishing, Cham (2019).https://doi.org/10.1007/978-3-030-31304-3 1[14] Murrugarra, D., Veliz-Cuba, A., Aguilar, B., Laubenbacher, R.: Identification of control targetsin boolean molecular network models via computational algebra. BMC Systems Biology (1),94 (2016). https://doi.org/10.1186/s12918-016-0332-x[15] Samaga, R., Kamp, A.V., Klamt, S.: Computing combinatorial intervention strategies andfailure modes in signaling networks. Journal of Computational Biology (1), 39–53 (2010).https://doi.org/10.1089/cmb.2009.0121[16] Takahashi, K., Yamanaka, S.: A decade of transcription factor-mediated reprogram-ming to pluripotency. Nature Reviews Molecular Cell Biology (3), 183–193 (2016).https://doi.org/10.1038/nrm.2016.8[17] Yang, G., G´omez Tejeda Za˜nudo, J., Albert, R.: Target control in logical mod-els using the domain of influence of nodes. Frontiers in Physiology , 454 (2018).https://doi.org/10.3389/fphys.2018.00454[18] Za˜nudo, J.G.T., Yang, G., Albert, R.: Structure-based control of complex networks with non-linear dynamics. Proceedings of the National Academy of Sciences (28), 7234–7239 (2017).https://doi.org/10.1073/pnas.1617387114[19] Za˜nudo, J.G.T., Albert, R.: Cell fate reprogramming by control of intracellular network dynamics.PLOS Computational Biology (4), 1–24 (2015). https://doi.org/10.1371/journal.pcbi.10041931320] Zhang, R., Shah, M.V., Yang, J., Nyland, S.B., Liu, X., Yun, J.K., Albert, R.,Loughran, T.P.: Network model of survival signaling in large granular lymphocyteleukemia. Proceedings of the National Academy of Sciences (42), 16308–16313 (2008).https://doi.org/10.1073/pnas.0806447105 Figure 5:
MAPK network, figure adapted from [7].14 able 1:
Apoptotic stable states of MAPK network (10). Each column represents a stable state. Thenumber in the cell indicates the value of the variable in the stable state. A A A A A A A A A A AKT 0 0 0 0 0 0 0 0 0 0AP1 1 1 1 1 1 1 1 1 1 1ATF2 1 1 1 1 1 1 1 1 1 1ATM 0 1 0 1 0 0 1 1 1 1Apoptosis 1 1 1 1 1 1 1 1 1 1BCL2 0 0 0 0 0 0 0 0 0 0CREB 1 1 1 1 1 1 1 1 1 1DNA-damage 0 1 0 1 0 0 1 1 1 1DUSP1 1 1 1 1 1 1 1 1 1 1EGFR 0 0 0 0 0 0 0 0 0 0EGFR-stimulus 1 1 0 0 1 0 1 0 0 0ELK1 1 1 1 1 1 1 1 1 1 1ERK 0 0 0 0 0 0 0 0 0 0FGFR3 0 0 0 0 0 0 0 0 0 0FGFR3-stimulus 1 1 1 1 0 0 0 0 0 0FOS 0 0 0 0 0 0 0 0 0 0FOXO3 1 1 1 1 1 1 1 1 1 1FRS2 0 0 0 0 0 0 0 0 0 0GAB1 1 1 1 1 1 1 1 1 1 0GADD45 1 1 1 1 1 1 1 1 1 1GRB2 1 1 1 1 1 1 1 1 0 0Growth-Arrest 1 1 1 1 1 1 1 1 1 1JNK 1 1 1 1 1 1 1 1 1 1JUN 1 1 1 1 1 1 1 1 1 1MAP3K1-3 1 1 1 1 1 1 1 1 0 0MAX1 1 1 1 1 1 1 1 1 1 1MDM2 0 0 0 0 0 0 0 0 0 0MEK1-2 0 0 0 0 0 0 0 0 0 0MSK 1 1 1 1 1 1 1 1 1 1MTK1 1 1 1 1 1 1 1 1 1 1MYC 1 1 1 1 1 1 1 1 1 1PDK1 1 1 1 1 1 1 1 1 1 0PI3K 1 1 1 1 1 1 1 1 1 0PKC 0 0 0 0 0 0 0 0 0 0PLCG 0 0 0 0 0 0 0 0 0 0PPP2CA 1 1 1 1 1 1 1 1 1 1PTEN 1 1 1 1 1 1 1 1 1 1Proliferation 0 0 0 0 0 0 0 0 0 0RAF 1 1 1 1 1 1 1 1 0 0RAS 1 1 1 1 1 1 1 1 0 0RSK 0 0 0 0 0 0 0 0 0 0SMAD 1 1 1 1 1 1 1 1 0 0SOS 1 1 1 1 1 1 1 1 0 0SPRY 0 0 0 0 0 0 0 0 0 0TAK1 1 1 1 1 1 1 1 1 0 0TAOK 0 1 0 1 0 0 1 1 1 1TGFBR 1 1 1 1 1 1 1 1 0 0TGFBR-stimulus 1 1 1 1 1 1 1 1 0 0p14 1 1 1 1 1 1 1 1 1 1p21 1 1 1 1 1 1 1 1 1 1p38 1 1 1 1 1 1 1 1 1 1p53 1 1 1 1 1 1 1 1 1 1p70 0 0 0 0 0 0 0 0 0 0 able 2: Control strategies up to size 5 obtained for the apoptotic stable states of the MAPK network.Each column represents a control strategy for the indicated stable state. A number in a cell indicatesthe value to which the variable is fixed in the control strategy. Empty cells denote uncontrolledcomponents. A A A A A A A A A DNA-damage 0 1 0 1 0 0 1 1 1 1EGFR-stimulus 1 1 0 0 1 0 1 0 0 0FGFR3-stimulus 1 1 1 1 0 0 0 0 0 0TGFBR-stimulus 1 1 1 1 1 1 1 1 0 0GAB1 1PI3K 1