Torsten Schaub

Anne Siegel

Santiago Videla

S Torsten

A Chaub

Nnesiegel

Santiagovidela

Reasoning on the

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Reasoning on the response of logical signaling networks with Answer Set Programming Logical networks provide a simple yet powerful qualitative modeling approach in systems biology. In what follows, we focus on modeling the response of logical signaling networks by means of automated reasoning using Answer Set Programming (ASP). In this context, the problem consisting of learning logical networks is crucial in order to achieve unbiased and robust discoveries. Furthermore, it has been shown that many networks can be compatible with a given set of experimental observations. Thus, first we discuss how ASP can be used to exhaustively enumerate all these logical networks. In practice, this is a key step for the design of new experiments in order to reduce the uncertainty provided by such a family of models. Next, in order to gain control over the system, we look for intervention strategies that force a set of target species into a desired steady state. Altogether, this constitutes a pipeline for reasoning on logical signaling networks providing robust insights to systems biologists. In this chapter, we illustrate the usage of ASP for solving the aforementioned problems and discuss the novelty of our approach with respect to existing methods.

Introduction

Systems biology and signaling networks

Systems biology is an interdisciplinary field aiming at the investigation and understanding of biology at a system and multi-scale level [IDE 01, KIT 02]. After biological entities have been identified in a specific environment, it remains to elucidate Chapter written by Torsten SCHAUB,AnneSIEGEL and Santiago VIDELA.

how they interact with each other in order to carry out a particular biological function. Thus, rather than focusing on the components themselves, one is interested in the nature of the links that connect them and the functionalities arising from such interactions. Notably, advances on high-throughput experimental technologies have been one of the main driving forces of systems biology. Such technologies have allowed biologists to study biological systems as a whole rather than as individual components. Nevertheless, the "reductionist" approach of molecular biology has been fundamental for the construction of the large cataloguesofbiologicalentitiesavailablenowadays. In fact, some authors have considered systems biology not as a new field of research, but instead as an approach to biomedical research combining "reductionist" and "integrationist" techniques [KOH 10].

As it is often the case, the application of novel technologies has led to profound conceptual and philosophical changes in biology. From the early days of molecular biology existed the idea that the DNA sequence dictates most of cell actions, as do the instructions in a computer program. Recently, together with the advent of systems biology, such a mechanistic understanding has been strongly revisited. Instead, an informatic perspective on the role of the genome has been established. From this point of view, the focus is on what the cell does with and to its genome products rather than on what the genome directs the cell to execute [SHA 09]. Then, for any biological system, one can envision at least a three-way interaction between DNA products, the environment and the phenotype [KOH 10]. In this scheme, the group of entities mediating between such interactions are the so-called biological networks. Deciphering the functioning of these complex networks is the central task of systems biology. Importantly, in order to cope with the increasing complexity of large-scale networks, mathematical and computational modeling is required. Hence, the development of such modeling approaches is a major goal in the field.

From the early millennium, many efforts have been made to develop relevant formalisms and modeling frameworks to take into account the specificities of complex biological systems. Among them, one can distinguish between mathematical and computational modeling approaches [FIS 07]. Essentially, mathematical (or quantitative) models are based on denotational semantics, that is, models are specified by mathematical equations describing how certain quantitites change over time. On the other hand, computational (or qualitative) models are based on operational semantics, that is, models are specified in terms of a sequence of steps describing how the states of an abstract machine relate to each other (not necessarily deterministically). Notably, each type of models provides a different levelofabstractionenablingtoaddressdifferent kinds of questions. Inf a c t ,h y b r i dm o d e l i n gp r e c i s e l ya i m sa te x p l o i t i n gt h e best of both worlds whenever possible. In any case, it is clear that intuition is not enough to face the complexity of large-scale biological systems. Thus, systematic and elaborated methodological tools are required by (systems) biologists. Moreover, the development of such modeling frameworks is leading to a hypothesis-driven research in biology [IDE 01, KIT 02]. At first, due to the lack of information, multiple hypotheses are usually generated from prior knowledge by either mathematical, computational or hybrid modeling. Next, decision-making methods can be used to suggest new experiments in order to reduce ambiguous hypotheses [KRE 09]. Finally, new experimental data is produced to test the generated hypotheses, models are refined, and the loop is started over again. Interestingly, to some extent, this iterative process could be automatized allowing an autonomous scientific discovery [SPA 10].

Among the biological networks mediating between genes, the environment, and the phenotype, signal transduction networks are crucial for the understanding of the response to external and internal perturbations. To be more precise, signal transduction occurs when an extracellular signaling molecule binds to a specific cell surface receptor protein. Such a binding causes a conformational change in the receptor that initiates a sequence of reactions leading to a specific cellular response such as growth, survival, apoptosis (cell death), and migration. Post-translational modifications, notably protein phosphorylation, play a key role in signaling. Importantly, signaling networks are involved in biomedical processes and their control has a crucial impact on drug target identification and diagnosis. Nowadays, there exist public repositories such as Pathways Commons [ CER 11], Pathways Interaction Database [SCH 09a], and KEGG [KAN 10] that contain curated knowledge about intracellular causal molecular interactions, from which canonical cell signaling networks can be retrieved [GUZ 12]. Such biological networks are derived fromv a s tg e n e r i ck n o w l e d g ec o m p i l e df r o m different cell types. Nevertheless, little is known about the exact chaining and composition of signaling events within these networks in specific cells and specific conditions. For example, in cancer cells, signaling networks frequently become compromised, leading to abnormal behaviors and responses to external stimuli. Many current and emerging cancer treatments are designed to block nodes in signaling networks, thereby altering signaling cascades. Thus, advancing our understanding of how these networks are deregulated across specific environments will ultimately lead to more effective treatment strategies for patients. In fact, there is emerging experimental evidence that the combinatorics of interactions between cellular components in signaling networks is a primary mechanism for generating highly specialized biological systems [PAP 05]. In this context, phosphorylation assays are a recent form of high-throughput data providing information about protein-activity modifications in as p e c i fi cc e l lt y p eu p o nv a r i o u sc o m b i n a t o r i a lp e r t u r b a t i o n s[ A L E1 0 ] . T h e r e f o r e , moving beyond the aforementioned causal and canonical interactions towards mechanistic and specialized descriptions of signaling networks is a major challenge in systems biology. Importantly, cellular signaling networks operate over a wide range of timescales (from fractions of a seconds to hours). Thus, taking this into account often leads to significant simplifications [PAP 05, MAC 12].

Finally, we conclude this brief introduction to systems biology and signaling networks by noting that biological functionality is multilevel [NOB 10]. That is, signaling networks by no means could function in isolation from metabolic and regulatory processes. Hence, integrative modeling approaches considering these multiple levels of causation pose indeed a long-term goal in the field.

Logical signaling networks: state of the art and current challenges

Among various discrete dynamical systems, Boolean networks [KAU 69] provide arelativelysimplemodelingapproachyetabletocaptureinterestingandrelevantbehaviors in the cell. Especially in the case of poorly understood biological systems where quantitative (fine-grained) information is often scarce and hard to obtain. In fact, numerous successful examples together with comprenhesive methodological reviews can be found in [RÉK 08, MOR 10, WAN 12, SAA 13]. In particular, it has been shown that the response in signaling networks can be appropriately modeled with Boolean networks, as illustrated on several signal transduction pathways involved in diverse processes such as proliferation, cellc y c l er e g u l a t i o n ,a p o p t o s i s ,o rd i f f e r e ntiation [SAE 07, SAM 09, SAE 11, CAL 10]. Nonetheless, a large majority of the authors working with (Boolean) logic-based models, rely on ad-hoc methods to build their models. That is, in most cases, Boolean networks are constructed manually based on information extracted from biological literature and available experimental data. Notably, the manual identification of logic rules underlying the system being studied is often hard, error-prone, and time consuming. In this context, researchers typically aim at modeling a given biological system by means of one Boolean network only. Afterwards, dynamical and structural analysis (often computationally demanding) are conducted over the model. Clearly, biological insights and novel hypotheses resulting from these analysis may be incomplete, incorrect or biased if the model was not precise enough. Therefore, automated learning of (Boolean) logic models is required in order to achieve unbiased and robust discoveries.

Reverse engineering in systems biology consists of building mathematical and computational models of biological systems based on the available knowledge and experimental data. Towards the construction of predictive models, one can convert the generic prior knowledge (for instance, canonical cell signaling networks) into a quantitative or qualitative model (for example, a set of differential equations or a set of logic rules) that can be simulated or executed. Next, if enough experimental data is available, the model can be fitted to it in order to obtain the most plausible models for certain environmental conditions or specific cell type. This is normally achieved by defining an objective fitness function to be optimized [BAN 08]. Optimization over quantitative modeling leads to continuous optimization problems. On the other hand, reverse engineering considering qualitative models typically give rise to combinatorial (discrete) optimization problems. Such inferences (or learning) of either quantitative or qualitative models, have been successfully applied for regulatory, signaling, and metabolic networks. Notably, it represents a very active area of research as illustrated by the successive "DREAM" challenges [STO 07].

In the context of logic-based models, the inference of Boolean (genetic) networks from time-series gene expression data has been addressed by several authors under different hypotheses and methods [LIA 98, AKU 00, IDE 00, LÄH 03]. Recently, a brief review and evaluation of these methods has been published in [BER 13]. Importantly, methods for reverse engineering of biological systems are highly dependent on available (amount of) data, prior knowledge and modeling hypotheses. In particular, reverse engineering of Boolean logic models by confronting prior knowledge on causal interactions with phosphorylation activities has been first described in [SAE 09]. A genetic algorithm implementation was proposed to solve the underlying optimization problem, and a software was provided, CellNOpt [TER 12]. Nonetheless, stochastic search methods cannot characterize the models precisely: they are intrinsically unable not just to provide a complete set of solutions, but also to guarantee that an optimal solution is found. To overcome some of this shortcomings, mathematical programming approaches were presented in [MIT 09, SHA 12]. Notably, authors in [SAE 09] have shown that the model is very likely to be non-identifiable when we consider the experimental error from measurements. Hence, rather than looking for the optimum Boolean model, one is interested in finding (nearly) optimal models within certain tolerance. Interestingly, in the context of quantitative modeling, authors in [CHE 09] have elaborated upon the same argument. Clearly, an exhaustive enumeration of (nearly) optimal solutions would allow for identifying admissible Boolean logic models without any methodological bias. Importantly, previous methods, namely stochastic search and mathematical programming, are not able to cope with this question. Moreover, all subsequent analysis will certainly profit from having such a complete characterization of feasible models. For example, finding "key-players" in signaling networks leading to potential therapeutic interventions is a highlyrelevantquestionforcancerresearchand biomedicine in general. In the context of logic-based models for signaling networks, this question has been recently addressed by [ABD 08, SAM 10, WAN 11,LAY 11]. Despite the specific problem settings and computational approaches in each work, all of them consider a single logic model describing the system. Moreover, due to computational limitations they are restricted to look for a small number (or even single) interventions. Notably, interventions allowing for accomplish certain goals in a given model are very likely to fail in another model which may describe the system similarly well. Therefore, being able to address the same question but considering a family of feasible models may lead to more robust solutions. In fact, this is in line with recent work showing that an ensemble of models often yields more robust predictions than each model in isolation [KUE 07,MAR 12]. In this context, there is an increasing demand of more powerful computational methods in order to achieve robust discoveries in systems biology.

Authors contribution

In this chapter, we present a generic, flexible, and unified framework for modeling logical networks and perform automated reasoning on their response. More precisely, we characterize the response of logical networks under either two-or three-valued logics based on fixpoint semantics. Furthermore, a representation using Answer Set Programming (ASP; [BAR 03]) is provided. ASP is a declarative problem solving paradigm, in which a problem is encoded as a logic program such that its so-called answer sets represent solutions to the problem. Furthermore, some of the key-aspects of ASP like its ease to express defaults and reachability are particularly relevant for dealing with biological networks. Importantly, available systems nowadays provide arichyetsimplemodelinglanguage,high-performance solving capacities, and automated reasoning modes [GEB 12a]. Altogether, ASP provides a powerful computational framework for addressing hard combinatorial problems in systems biology by reasoning over the complete search space. Notably, our ASP representation can be easily elaborated in order to consider specific problem settings as we illustrate in subsequent sections. In particular, in the context of signaling networks, we address the problem consisting of learning Boolean logic models of immediate-early response and the problem consisting of finding minimal intervention strategies in logical networks.

The remainder of the chapter is structured as follows: Section 1.2 introduces ASP and its main features; Section 1.3 describes our framework for modeling logical networks and two specific problems with their corresponding modeling and solving using ASP systems; and Section 1.4 concludes and poses prospective challenges.

Answer Set Programming at a glance

Answer Set Programming (ASP; [BAR 03, GEB 12a]) provides a declarative framework for modeling combinatorial problems in Knowledge Representation and Reasoning. The unique pairing of declarativeness and performance in state-of-the-art ASP solvers allows for concentrating on an actual problem, rather than a smart way of implementing it. The basic idea of ASP ist oe x p r e s sap r o b l e mi nal o g i c a lf o r m a t so that the models of its representation provide the solutions to the original problem. Problems are expressed as logic programs and the resulting models are referred to as answer sets. Although determining whether a program has an answer set is the fundamental decision problem in ASP, more reasoning modes are needed for covering the variety of reasoning problems encountered in applications. Hence, a modern ASP solver, like clasp [GEB 12b] supports several reasoning modes for assessing the multitude of answer sets, among them, regular andp r o j e c t i v ee n u m e r a t i o n ,i n t e r s e c t i o n and union, and multi-criteria optimization. As well, these reasoning modes can be combined, for instance, for computing the intersection of all optimal models. This is accomplished in several steps. At first, a logicp r o g r a mw i t hfi r s t -o r d e rv a r i a b l e si s turned into a propositional logic program by means of efficient database techniques. This is in turn passed to a solver computing the answer sets of the resulting program by using advanced Boolean constraint technology. For optimization, a solver like clasp uses usually branch-and-bound algorithms but other choices, like computing unsatisfiable cores, are provided as well. The enumeration of all optimal models is done via the option --opt-mode=optN.Atfirstanoptimalmodelisdeterminedalongwith its optimum value. This computation has itself two distinct phases. First, an optimal model candidate must be found and second, it must be shown that there is no better candidate; the latter amounts to a proof of unsatisfiability and is often rather demanding (because of its exhaustive nature). Then,allmodelspossessingtheoptimumscore are enumerated. Notice that this way one can enumerate all (strictly) optimal solutions. Nonetheless, we are often interested in (nearly) optimal answer sets as well. For a concrete example on how we address this in practice, we refer the reader to the encoding provided in Listing 1.8 and its solving in Listing 1.10.

Our encodings are written in the input language of gringo 4s e r i e s . S u c hal a nguage implements most of the so-called ASP-Core-2 standard. 1 In what follows, we introduce its basic syntax and we refer the reader to the available documentation for more details. An atom is a predicate symbol followed by a sequence of terms (e.g. p(a,b),q(X,f(a,b))). A term is a constant (e.g. c, 42)o raf u n c t i o ns y m b o lf o llowed by a sequence of terms (e.g. f(a,b), g(X,10))whereuppercaselettersdenote first-order variables. Then, a rule is of the form

h:-b 1 ,...,b n .
where h (head) is an atom and any b j (body) is a literal of the form a or not a for an atom a where the connective not corresponds to default negation. The connectives :-and , can be read as if and and,r e s p e c t i v e l y . F u r t h e r m o r e ,ar u l ew i t h o u tb o d y is a fact,w h e r e a sar u l ew i t h o u th e a di sa nintegrity constraint.A l o g i c p r o g r a m consists of a set of rules, each of which is terminated by a period. An atom preceded with default negation, not,i ss a t i s fi e du n l e s st h ea t o mi sf o u n dt ob et r u e . I nA S P , the semantics of a logic program is given by the stable models semantics [GEL 88]. Intuitively, the head of a rule has to be true whenever all its body literals are true. This semantics requires that each true atom must also have some derivation, that is, an atom cannot be true if there is no rule deriving it. This implies that only atoms appearing in some head can appear in answer sets.

We end this quick introduction by three language constructs particularly interesting for our encodings. First, the so called choice rule of the form,

{h 1 ; ...; h m }:-b 1 ,...,b n .
allows us to express choices over subsets of atoms. Any subset of its head atoms can be included in an answer set, provided the body literals are satisfied. Note that using a choice rule one can easily generate an exponential search space of candidate solutions. Second, a conditional literal is of the form l : l 1 ,...,l n

The purpose of this language construct is to govern the instantiation of the literal l through the literals l 1 ,...,l n .I n t h i s r e s p e c t ,t h e c o n d i t i o n a l l i t e r a l a b o v e c a n b e regarded as the list of elements in the set {l | l 1 ,...,l n }.F i n a l l y ,f o rs o l v i n g( m u l t icriteria) optimization problems, ASP allows for expressing (multiple) cost functions in terms of a weighted sum of elements subject to minimization and/or maximization. Such objective functions are expressed in gringo 4intermsof(se veral)optimization statements of the form

#opt{w 1 @l 1 ,t 11 ,...,t m1 : b 1 ,...,b n1 ; ...; w k @l k ,t 1 k ,...,t m k : b 1 k ,...,b n k }.
where #opt ∈ {"#minimize", "#maximize"}, w i ,l i ,t 1i ,...,t mi are terms and b 1i ,...,b ni are literals for k ≥ 0, 1 ≤ i ≤ k, m i ≥ 0 and n i ≥ 0.F u r t h e r m o r e , w i and l i stand for an integer weight and priority level.P r i o r i t i e sa l l o wf o rr e p r esenting lexicographically ordered optimization objectives, greater levels being more significant than smaller ones.

Answer Set Programming for Systems biology

Our work contributes to a growing list of ASP applications in systems biology. Almost a decade ago, Baral et Let f : X → Y ,bea(partial)functionmappingvaluesx ∈ X ′ ⊆ X to values y ∈ Y .W ede not et hes e tofv a l ue sx such that f (x) is defined, i.e. X ′ ,withdom(f ).W e sometimes represent mappings extensionally as sets, viz. {x → f (x) | x ∈ dom(f )}, for checking containment, difference, etc. To avoid conflicts when composing truth assignments, we define

A • B =(A \ B) ∪ B where B = {v → s | v → s ∈ B} and t = f , f = t, u = u.

Reasoning on the response of logical networks

Logical networks

A logical network consists of a finite set V of propositional variables and a (partial) function φ mapping a variable v ∈ V to a propositional formula φ(v) over V . The logical steady states of (V, φ) are given by truth assignments yielding identical values for v and φ(v) for all v ∈ dom(φ).G e n e r a l l ys p e a k i n g ,s u c hl o g i c a ln e t w o r k s can be seen as synchronous Boolean networks [KAU 69]. However, since we consider both, two-and three-valued logics, we refrain from using the term "Boolean". Without loss of generality, we assume only formulas in disjunctive normal form.2 For illustration, let us consider the logical network consisting of the set V of species variables {i 1 ,i 2 ,a,...,g,o 1 ,o 2 } along with the function φ defined as:

a ∧ b b t f u a t t f u f f f f u u f u a ∨ b b t f u a t t t t f t f u u t u u a ¬a t f f t u u
φ = ⎧ ⎨ ⎩ a → ¬d b → a ∧ i 1 c → b ∨ e d → c e → ¬i 1 ∧ i 2 f → e ∨ g g → f o 1 → c o 2 → g ⎫ ⎬
⎭ Note that φ leaves the specification of the variables i 1 and i 2 undefined. Furthermore, we represent logical networks as (signed) directed hypergraphs as shown in Figure 1.1(a). A (signed) directed hypergraph is defined by a pair (V, H) with vertices V and (signed) directed hyperedges H;a(s i gned)di rect edhyperedgei sapai r(S, t) where S is a finite, non-empty set of pairs 3 Then, we say that the (signed) directed hypergraph (V, H) represents the logical network (V, φ) if and only if for every v ∈ dom(φ) and variable w ∈ V that occurs positively (resp. negatively) in some conjunct ψ of φ(v),thereisahyperedge (S ψ ,v) with (w, 1) ∈ S ψ (resp. (w, -1) ∈ S ψ ); and vice versa. Following the example shown in Figure 1.1(a), if we consider the mapping φ(e)=¬i 1 ∧ i 2 ,weneedto verify the existence of the hyperedge (S ¬i1∧i2 ,e) with S ¬i1∧i2 = {(i 1 , -1), (i 2 , 1)}. Similarly, for the mapping φ(c)=b ∨ e,weneedtoverifytheexistenceofthehyperedges (S b ,c) with S b = {(b, 1)} and (S e ,c) with S e = {(e, 1)}.

(v i ,s i ) with v i ∈ V, s i ∈ {1, -1} and t ∈ V .
Next, we introduce the notion of the interaction graph underlying a logical network

(V, φ).A ninteraction graph (V, E, σ) is a signed and directed graph with vertices V , directed edges E ⊆ V × V and signature σ ⊆ E × {1, -1}.M o r e o v e r ,w es a yt h a t Σ (V,φ) =( V, E, σ) is the underlying interaction graph of (V, φ) if for every edge (v, w) ∈ E with ((v, w), 1) ∈ σ (resp. ((v, w), -1) ∈ σ)
, the variable v occurs positively (resp. negatively) in the formula φ(w).N o t et h a tt h e r ei saone-to-many relation in the sense that the same graph (V, E, σ) corresponds to the underlying interaction graph Σ (V,φ) for possibly many logical networks (V, φ).N o w ,w ec a nr e l yo ns t a n d a r d notions from graph theory to capture severalconceptsonlogicalnetworks.Recallthat apathinagraphisasequenceofedgesconnectingasequenceofvertices.Thelength of a path is given by the number of edges whereas its sign is the product of the signs of the traversed edges. Herein, we consider only paths with length greater than zero. Thus, an edge (v, v) is required in order to consider the existence of a path from v to v.W es a yt h e r ei sap o s i t i v e( r e s p .n e g a t i v e )p a t hf r o mv to w in (V, φ) if and only if there is a positive (resp. negative) path from v to w in Σ (V,φ) .F u r t h e r m o r e ,w e say there is a positive (resp. negative) feedback-loop in (V, φ) if and only if for some v ∈ V there is a positive (resp. negative) path from v to v in Σ (V,φ) .

Characterizing the response of the system

Let (V, φ) be a logical network describing a biological system of interest. For capturing changes in the environment of such a biological system, for instance, due to an experimental intervention (over-expression or knock-out), we introduce the notion of clamping variables in the network for overriding their original specification. To this end, we define a clamping assignment C as a partial two-valued truth assignment over V .T hen,w edefi net hemappi ngφ| C as

φ| C (v)= ⎧ ⎨ ⎩ ⊤ if v → t ∈ C ⊥ if v → f ∈ C φ(v) if v ∈ dom(φ) \ dom(C) yielding the modified logical network (V, φ| C ).M o r e o v e

r , i t i s w o r t h n o t i n g t h a t dom(φ) ⊆ dom(φ| C ).L e tusi l l us t r a t et hi sw i t hourt o ye xa mpl ei nF i gur e1. 1( a ) .L e t C be the clamping assignment defined by {i

1 → t,i 2 → f ,g → f }.T h e n ,φ| C is a complete mapping over V defined as: φ| C = ⎧ ⎨ ⎩ i 1 →⊤ i 2 →⊥ a → ¬d b → a ∧ i 1 c → b ∨ e d → c e → ¬i 1 ∧ i 2 f → e ∨ g g →⊥ o 1 → c o 2 → g ⎫ ⎬ ⎭
In practice, clamping assignments are usually restricted to a subset of variables X ⊆ V .M o r e o v e r ,c e r t a i nv a r i a b l e si nX may be further restricted to be clamped either to a single truth value, viz. t or f ,o rn o tc l a m p e da ta l l . T h e s er e s t r i c t i o nw i l l be typically related to context-specific application settings, for instance, the kind of biological entity described by each variable and "real-world" experimental limitations over such entity.

Next, for capturing the synchronous updates in a logical network (V, φ),w ef o llow [ INO 11] and define the single-step operator on either two-or three-valued (complete) truth assignments over V :4

Ω (V,φ) (A)={v → A(φ(v)) | v ∈ dom(φ)} ∪ {v → A(v) | v ∈ V \ dom(φ)} ,
where A is extended to formulas in the standard way. Notice that the definition above captures the fact that unmapped variables in φ remain unchanged with respect to the value assigned in A.F u r t h e r m o r e ,f o rc a p t u r i n gt h et r a j e c t o r yo fs t a t eA we define the iterative variant of Ω (V,φ) as

Ω 0 (V,φ) (A)=A and Ω j+1 (V,φ) (A)=Ω (V,φ) Ω j (V,φ) (A) .
In biological terms, a sequence (Ω j (V,φ) (A)) j∈J captures the signal propagation starting in state A.I np a r t i c u l a r ,w ea r ei n t e r e s t e di nt h efi x p o i n to fΩ (V,φ) reachable from certain initial assignment A.I m p o r t a n t l y ,t h ee x i s t e n c eo fs u c hafi x p o i n ti sn o tn e cessarily guaranteed. In general, it depends on the definition of A and the presence or absence of feedback-loops in (V, φ).B u ti nc a s eo fe x i s t e n c e ,s u c hafi x p o i n td escribes a logical steady state which is interpreted as the response of the biological system described by (V, φ).T ob em o r ep r e c i s e ,t h ec h o i c eo fA is related to how we model the absence of information in the context of either two-or three-valued logics. Hence, when we consider three-valued logic, we use the initial assignment

A u = {v → u | v ∈ V }.I n t e
r e s t i n g l y ,i nt h i sc o n t e x tafi x p o i n ti sr e a c h e dr e g a r dless of the presence or absence of feedback-loops in the network. Moreover, such a fixpoint poses the property that each of its variables is assigned to u unless there is a cause to assign it to either t or f .O nt h eo t h e rh a n d ,w h e nw ec o n s i d e rt w o -v a l u e d logic, we use the initial assignment

A f = {v → f | v ∈ V }.U n f o r t u
n a t e l y ,i nt h i s context, the presence of feedback-loops typically avoids reaching a fixpoint. Next, let us illustrate the iterated application of Ω (V,φ|C) for our toy example in the context of both, two-and three-valued logics. Recall that we have defined φ| C above for clamp-

ing assignment C = {i 1 → t,i 2 → f ,g → f }.T h er e s u l t i n ga s s i g n m e n t sf r o mt h e computation of Ω j (V,φ|C) (A) with either A = A u or A = A f are shown in Table 1.3. Notably, when we consider three-valued logic, Ω 3 (V,φ|C) (A u )=Ω 4 (V,φ|C) (A u ) results
in the fixpoint:

⎧ ⎨ ⎩ i 1 → t i 2 → f a → u b → u c → u d → u e → f f → f g → f o 1 → u o 2 → f ⎫ ⎬ ⎭ Meanwhile, under two-valued logic we obtain Ω 1 (V,φ|C) (A f )=Ω 9 (V,φ|C ) (A f
) which leads to an oscillatory behavior for variables a, b, c, d and o 1 .N o t i c et h a tt h e s ev a r iables correspond to the ones assigned to u in the fixpoint reached for Ω 3 (V,φ|C) (A u ). In this case, the oscillatory behavior is induced by the negative feedback-loop over the path a, b, c, d.T h u s ,o n ec a nv e r i f yt h a t ,f o re x a m p l e ,i fw er e m o v et h em a p p i n g d → c from the definition of φ (leaving d undefined in φ), then Ω 4 (V,φ|C) (A f ) would result in the fixpoint:

⎧ ⎨ ⎩ i 1 → t i 2 → f a → t b → t c → t d → f e → f f → f g → f o 1 → t o 2 → f ⎫ ⎬ ⎭
In fact, whenever we consider a logical network (V, φ) under two-valued logic, we enforce that (V, φ) is free of feedback-loops. Notably as detailed below, although not capable of capturing dynamical properties, this simplification guarantees the existence of a fixpoint while it allows us to characterize the so-called iimmediate-early response in signaling networks.

Logical networks and their response with Answer Set Programming

Let (V, φ) be a logical network. We represent the variables V as facts over the predicate variable/1,namely,variable(v) for all v ∈ V .R e c a l lt h a tw ea s s u m e φ(v) to be in disjunctive normal form for all v ∈ V .H e n c e ,φ(v) is a set of clauses and a clause a set of literals. We represent formulas using predicates formula/2, dnf/2,a n dclause/3.T h ef a c t sformula(v,s φ(v) ) map variables v ∈ V to their corresponding formulas φ(v),f a c t sdnf(s φ(v) ,s ψ ) associate φ(v) with its clauses ψ ∈ φ(v),f a c t sclause(s ψ ,v,1) associate clause ψ with its positive literals v ∈ ψ ∩ V ,andfactsclause(s ψ ,v,-1) associate clause ψ with its negative literals ¬v ∈ ψ.N o t e t h a t e a c hs (•) stands for some arbitrary but unique name in its respective context. Listing 1.1 shows the representation of our toy example logical network in Figure 1.1(a).

Listing 1.1: Logical networks representation as logical facts variable (i1). variable (i2). variable (o2). variable (o1). variable (a). variable (b). variable (c). variable (d). variable (e). variable (f). variable (g).

formula (a,0). formula(b,2). formula (c,1). formula (d,4). formula (e,3). formula(f,6). formula (g,5). formula (o1 ,4). formula (o2 ,7). dnf (0,5). dnf(1,6). dnf (1,0). dnf (2,3). dnf(3,7). dnf (4,1). dnf(5,2). dnf (6,4). dnf (6,6). dnf(7,4). clause(0,b,1). clause (1,c,1). clause(2,f,1). clause(3,a,1). clause(3,i1,1). clause (4,g,1). clause(5,d,-1). clause(6,e,1). clause(7,i2,1). clause (7,i1,-1).

The representation of clamping assignments is straightforward. Note that in the following we use 1 and -1 for truth assignments to t and f ,r e s p e c t i v e l y . L e tC be ac l a m p i n ga s s i g n m e n to v e rV ,w er e p r e s e n tt h ea s s i g n m e n t si nC as facts over the predicate clamped/2,namely ,clamped(v,s) with

s =1if C(v)=t and s = -1 if C(v)=f .T h ee x a m p l ec l a m p i n ga s s i g n m e n tC = {i 1 → t,i 2 → f ,g → f } is shown in Listing 1.2.
Listing 1.2: Clamping assignment as logical facts clamped (i1 ,1). clamped (i2,-1). clamped (g,-1).

Furthermore, we introduce two rules deriving predicates eval/2 and free/2.T h e predicate eval/2 captures the fact that clamped variables are effectively fixed to the Next, we describe how we model either two-or three-valued logics in ASP. In fact, the rule modeling the propagation of (positive) true values is the same for both logics. Essentially, we exploit the fact that formulas φ(v) are in disjunctive normal form. Hence, under both logics we derive eval(v,1) if v is not clamped and there exists a conjunct ψ ∈ φ(v) such that all its literals evaluate positively. The rule describing this is shown in Listing 1.4.

φ| C (v) ⊤⊥¬da∧ i 1 b ∨ ec¬i 1 ∧ i 2 e ∨ g ⊥ cg v ∈ V i 1 i 2 ab c d e f g o 1 o 2 3-valued Ω 0 (V,φ|C ) (A u ) uu u u u u u u uu u Ω 1 (V,φ|C ) (A u ) tf u u u u u u fuu Ω 2 (V,φ|C ) (A u ) tf u u u u f u fuf Ω 3 (V,φ|C ) (A u ) tf u u u u f f fuf Ω 4 (V,φ|C ) (A u ) tf u u u u f f fuf 2-valued Ω 0 (V,φ|C) (A f ) ff f f f f f f ff f Ω 1 (V,φ|C) (A f ) tf t f f f f f ff f Ω 2 (V,φ|C) (A f ) tf t t f f f f ff f Ω 3 (V,φ|C) (A f ) tf t t t f f f ff f Ω 4 (V,φ|C) (A f ) tf t t t t f f f t f Ω 5 (V,φ|C) (A f ) tf f t t t f f f t f Ω 6 (V,φ|C) (A f ) tf f f t t f f f t f Ω 7 (V,φ|C) (A f ) tf f f f t f f f t f Ω 8 (V,φ|C) (A f ) tf f f f f f f ff f Ω 9 (V,φ|C) (A f ) tf t f f f f f ff f Table 1.3. Exemplary iterated application of Ω (V,φ| C ) for (V, φ) in Figure 1.1(a), clamping assignment C = {i1 → t,i2 → f ,c → f }
Listing 1.4: Positive propagation common to two-and three-valued logics eval(V,1) :-free(V,D); eval(W,T) : clause(J,W,T); dnf(D,J).

Meanwhile, the propagation of (negative) false values depends on the type of logic under consideration. On the one hand, when we consider two-valued logic, we use the rule shown in Listing 1.5 to derive eval(v,-1) if it cannot be proved that v evaluates positively, that is, not eval(v,1).

Listing 1.5: Negative propagation for two-valued logic (with default negation) eval(V,-1) :-variable (V); not eval(V,1).

On the other hand, when we consider three-valued logic, we use the rules shown in Listing 1.6. Notice that in this case, we derive eval(v,-1) only if it can be proved that all clauses ψ ∈ φ(v) evaluate negatively. A clause φ evaluates negatively if at least one of its literals evaluates negatively. Clauses evaluating negatively are represented with the predicate eval_clause/2. Listing 1.6: Negative propagation for three-valued logic (with explicit proof) eval_clause (J,-1) :-clause(J,V,S); eval(V,-S). eval(V,-1) :-free(V,D); eval_clause (J,-1) : dnf(D,J).

Interestingly, our ASP representation is relatively simple yet flexible enough to be extended and adapted for specific applications as we illustrate in the remainder of this chapter. In what follows, we provide a formal characterization for two very relevant problems over logical networks together with the corresponding modeling and solving based on the presented ASP encodings. First we extend our representation to learn logical networks from a given interaction graph confronting their response with experimental observations. Moreover, we consider several clamping assignments (describing experimental conditions) simultaneously instead of only one. Afterwards, we adapt our representation again aiming at reasoning over a family of logical networks and finding clamping assignments (describing therapeutic interventions) leading to responses satisfying specific goals.

Learning models of immediate-early response

Background

Firstly, we briefly summarize the main biological hypotheses in [SAE 09] providing the foundation for the concept of Boolean logic models of immediate-early response. Concretely, a Boolean logic model of immediate-early response is a logical network (V, φ) as defined above, without feedback-loops and using classical (Boolean) logics. The main assumption under Boolean logic models of immediateearly response is the following. The response of a biological system to external perturbations occurs at several time scales [PAP 05]. Thus, one can discriminate between fast and slow events. Under this assumption, at a given time after perturbation, the system reaches a state on which fast events are relevant, but slow events (such as protein degradation) have a relatively insignificant effect. In this context, we say that the system has reached a pseudo-steady state describing the early events or immediate-early response. Qualitatively, these states can be computed as logical steady states in the Boolean network (V, φ) [KLA 06b]. In fact, the discrimination between fast and slow events has an important consequence. Since we focus on fast or early events, it is assumed that oscillation or multi-stability caused by feedbackloops [REM 08, PAU 12] cannot happen until the second phase of signal propagation occurring at a slower time scale. Therefore, eedback-loops are not included in Boolean logic models of immediate-early response assuming that they will become active in a late phase [MAC 12]. Notably, it follows that starting from any initial state, aB oole a nlogicmode lofi mme dia te -e a r lyresponse reaches a unique steady state or fixpoint in polynomial time [PAU 12]. Thus, such modeling approach, although not capable of capturing dynamical properties, provides a relatively simple framework for input-output predictive models.

Based on the assumptions and concepts described above, authors in [SAE 09] have proposed a method to learn from an interaction graph and phosphorylation activities at a pseudo-steady state, Boolean logic models of immediate-early response fitting experimental data. In the remaining of this section we provide a precise characterization of this problem using the notions introduced in Section 1.3 and adapting our Answer Set Programming representation accordingly.

Problem

A prior knowledge network is an interaction graph (V, E, σ) as defined above. In addition, we distinguish three special subsets of species in V namely, the stimuli (V S ), the knock-outs (V K )andthereadouts (V R ). Nodes in V S denote extracellular ligands that can be stimulated or knocked-in and thus, we assume they have indegree equal to zero. Nodes in V K denote intracellular species that can be inhibited or knocked-out by various experimental tools such as small-molecule drugs, antibodies, or RNAi. Finally, nodes in V R denote species that can be measured by using an antibody. Notably, species in none of these sets, are neither measured, nor manipulated for the given experimental setup. Let us denote with V U the set of such nodes. Then, except for V R and V K that may intersect, the sets V S ,V K ,V R and V U are pairwise mutually disjoint. An early simplification consists on compressing the PKN in order to collapse most of the nodes in V U .T h i so f t e nr e s u l t so nas i g n i fi c a n tr e d u c t i o no ft h es e a r c hs p a c et h a t must be explored during learning. Thus, herein we assume a compressed PKN as an input and we refer the interested reader to [SAE 09] for a detailed description on this subject.

Given a PKN (V, E, σ),t h ec o n c e p to fa nexperimental condition over (V, E, σ) is captured by a clamping assignment over variables V S ∪ V K .R e c a l lt h a tc l a m ping assignments were defined above as partial two-valued assignments. To be more precise, while variables in V S can be clamped to either t or f ,v a r i a b l e si nV K can only be clamped to f .N e x t ,i fC is an experimental condition and v ∈ V S ,t h e n C(v)=t (resp. f )i n d i c a t e st h a tt h es t i m u l u sv is present (resp. absent), while if v ∈ V K ,t h e nC(v)=f indicates that the species v is inhibited or knocked out. In fact, since extracellular ligands by default are assumed to be absent, for the sake of simplicity we can omit clampings to f over variables in V S .T h e r e f o r e ,i fC is an experimental condition and v ∈ dom(C) then, either v ∈ V S and C(v)=t,o r v ∈ V K and C(v)=f .F u r t h e r m o r e ,t h ec o n c e p to fa nexperimental observation under an experimental condition C is captured by a partial mapping

P C : V R → [0, 1].
That is, dom(P C ) ⊆ V R denotes the set of measured readouts under the experimental condition C.I fv ∈ dom(P C ),t h e nP C (v) represents the phosphorylation activity at a pseudo-steady state of the readout v under C.N o t a b l y , i t i s r a t h e r c r i t i c a l t o choose a time point that is characteristic for the fast or early events in the biological system under consideration [MAC 12]. Since phosphorylation assays represents an average across a population of cells, the phosphorylation activity for each readout is usually normalized to [0, 1].F i n a l l y ,a nexperimental dataset ξ is a finite set of pairs (C i ,P Ci ) with experimental conditions C i and experimental observations P Ci .F u rther, we denote with N ξ the size of ξ given by the number of measured readouts across all experimental conditions i =1,...,n,i.e.,N ξ = n i=1 |dom(P Ci )|. Let us illustrate the concepts described above with our toy example. Consider the PKN (V, E, σ) defined in Figure 1.2(a). From the graph coloring, we have V S = {a, b, c}, V K = {d} and V R = {f, g}.Furthermore,letξ =((C 1 ,P C1 ),...,(C 4 ,P C4 )) be an example experimental dataset over (V, E, σ) defined by

C 1 = {a → t,c → t} P C1 = {f → 0.9,g → 0.0} C 2 = {a → t,c → t,d → f } P C2 = {f → 0.1,g → 0.9} C 3 = {a → t} P C3 = {f → 0.0,g → 0.1} C 4 = {a → t,b → t} P C4 = {f → 1.0,g → 0.8}.
(1.1)

In words, the experimental conditions C 1 ,...,C 4 can be read as follows. In C 1 ,stimuli a and c are present, stimulus b is absent and d is not inhibited; in C 2 ,stimulia, b, c are like in C 1 but d is inhibited; in C 3 ,o n l yt h es t i m u l u sa is present and d is not inhibited; and in C 4 ,stimulia and b are present, stimulus c is absent and d is not inhibited. Experimental observations P C1 ,...,P C4 give (normalized) phosphorylation activities for readouts f and g under the corresponding experimental condition.

Next, we introduce the notion of Boolean predictions. Let (V, E, σ) be a PKN. Further, let ξ =(C i ,P Ci ) be an experimental dataset over (V, E, σ) with i =1,...,n.

As detailed above, a Boolean logic model of immediate-early response is defined by al o g i c a ln e t w o r k(V, φ) without feedback-loops and using classical (Boolean) logics. Hence, now we can define the predictions (output) provided by a Boolean logic model of immediate-early response with respect to a given set of experimental conditions (input). Towards this end, we characterize the response of logical networks using fixpoint semantics as detailed above. More precisely, for i =1 ,...,n let F i be the fixpoint of Ω (V,φ|C i ) reachable from A f .N o t i c et h a ts u c hafi x p o i n ta l w a y s exists given that (V, φ) is free of feedback-loops. In words, each F i describe the logical response (starting from A f )of(V, φ) with respect to the experimental condition or clamping assignment C i .N e x t ,w ed e fi n eas t r a i g h t f o r w a r dt r a n s f o r m a t i o nf r o m truth values to binary but numerical values. Such a transformation provides a more convenient notation in order to compare predictions and phosphorylation activities. The Boolean prediction of (V, φ) with respect to the experimental condition C i is a function π i : V → {0, 1} defined as

π i (v)= 1 if F i (v)=t 0 if F i (v)=f .
As we see below, during learning we aim at explaining the given experimental dataset ξ.T h e r e f o r e ,w ea r ep a r t i c u l a r l yi n t e r e s t e do np r e d i c t i o n sw i t hr e s p e c tt ot h ee x p e r imental conditions included in ξ and over the measured variables in each experimental condition. Nevertheless, predictions with respect to non-performed experimental conditions and/or over non-observed species can be useful to generate testeable hypotheses on the response of the system.

As an example, consider the Boolean logic model (V, φ) from Figure 1.2(c) and the experimental condition C 2 from the dataset given in (1.1). Then, the clamped logical network (V, φ| C2 ) is defined by the mapping

φ| C2 = {a →⊤; b →⊥; c →⊤; d →⊥; e → b ∨ c; f → d ∧ e; g → e ∧ ¬c}.
Next, the fixpoint of Ω (V,φ|C 2 ) reachable from A f can be computed yielding the assignment F 2 defined as

F 2 = {a → t,b → f ,c → t,d → f ,e → t,f → f ,g → f }.
Finally, the Boolean prediction for (V, φ) under the experimental condition C 2 is given by

π 2 = {a → 1,b → 0,c → 1,d → 0,e → 1,f → 0,g → 0}.
We aim at learning Boolean logic models from a PKN and an experimental dataset. In fact, any learned model has to be supported by some evidence in the prior knowledge. To be more precise, given a PKN (V, E, σ) we consider only Boolean logic models (V, φ) without feedback-loops and such that, for each variable v ∈ V ,i fw occurs positively (resp. negatively) in φ(v) then, there exists an edge (w, v) ∈ E and ((w, v), 1) ∈ σ (resp. ((w, v), -1) ∈ σ). Towards this end, we consider a pre-processing step where the given PKN is expanded to generate a (signed) directed hypergraph describing all plausible logical interactions. For each v ∈ V having nonzero indegree, let Pred(v) be the set of its (signed) predecessors, namely, Pred(v)= {(u, s) | (u, v) ∈ E, ((u, v),s) ∈ σ}.F u r t h e r m o r e , l e t P(v) be the powerset of Pred(v),n a m e l y ,2 Pred(v) .T h e n , (V, H) is the (signed) directed hypergraph expanded from (V, E, σ) with nodes V and (signed) directed hyperedges H if for each v ∈ V , (p, v) ∈ H whenever p ∈ P(v).N e x t ,B o o l e a nl o g i cm o d e l sm u s te s s e n t i a l l y result from pruning (V, H).A d d i t i o n a l l y ,w ei m p o s et w oc o n s t r a i n t sr e l a t e dt ot h e fact that our Boolean logic models essentially aim at providing a framework for inputoutput predictions. 5 Firstly, for any variable u defined in φ all variables w ∈ φ(u) must be reachable from some stimuli variable. That is, we consider only Boolean logic models (V, φ) such that for every u ∈ dom(φ) and w ∈ φ(u),eitherw ∈ V S or there exist v ∈ V S and a path from v to w in the underlying interaction graph Σ (V,φ) . Secondly, every variable u defined in φ must reach some readout variable. That is, we consider only Boolean logic models (V, φ) such that for every u ∈ dom(φ) there exist v ∈ V R and a path from u to v in the underlying interaction graph Σ (V,φ) .F i n a l l y , let us denote with M (V,E,σ) the search space of Boolean logic models satisfying the conditions given above: evidence in (V, E, σ),n of e e d b a c k -l o o p s ,a n dr e a c h a b i l i t y from/to stimuli/readouts.

In Figure 1.2(a) we show an exemplary PKN and the corresponding expanded (signed) directed hypergraph in Figure 1.2(b). As already described, (signed) directed hypergraphs can be directly linked to logical networks. Thus, by considering each (signed) directed hyperedge in Figure 1.2(b) as either present or absent (and verifying the additional constraints related to feedback-loops and reachability from/to stimuli/readouts), one can generate the search space of Boolean logic models M (V,E,σ) defined above.

For a given PKN (V, E, σ),thereareexponentiallymanycandidateBooleanlogic models (V, φ) having an evidence on it. Therefore, authors in [SAE 09] put forward the idea of training Boolean logic models by confronting their corresponding Boolean predictions with phosphorylation activities atap s e u d o -s t e a d ys t a t e . I nt h i sc o n t e x t , two natural optimization criteria arise in order to conduct the learning: (1) model accuracy (biologically meaningful), and (2) model complexity (Occam's razor principle). In fact, this is a typical scenario on automatized learning of predictive models [FRE 04].

5. The interested reader may notice the analogy with the elimination of nodes neither controllable nor (leading to) observable during the compression of the PKN described in [SAE 09].

We now provide the precise formulationforeachoptimizationcriteria. Let(V, E, σ) be a PKN. Let ξ =( C i ,P Ci ) be an experimental dataset over (V, E, σ) with i = 1,...,n.L e t(V, φ) be a Boolean logic model having evidence in (V, E, σ) and let π 1 ,...,π n be its Boolean predictions with each π i defined under C i .Firstly ,basedon the residual sum of squares (RSS) we define the residual (Θ rss )of(V, φ) with respect to ξ as

Θ rss ((V, φ), ξ)= n i=1 v∈dom(PC i ) (P Ci (v) -π i (v)) 2 .
(1.2)

Secondly, for a given logical formula φ(v),l e tu sd e n o t ei t sl e n g t hb y|φ(v)|.T h e n , we define the size (Θ size )of(V, φ) as

Θ size ((V, φ)) = v∈dom(φ) |φ(v)|. (1.3)
Ap o p u l a ra n dr e l a t i v e l ys i m p l ea p p r o a c ht ocope with multi-objective optimization is to transform it into a single-objective optimization. Towards this end, one usually combines all criteria by defining a function using free parameters in order to assign different weights to each criteria. In fact, this is exactly the approach adopted in [SAE 09]. Therein, a single-objective function is defined that balances residual and size using a parameter α chosen to maximize the predictive power of the model. Moreover, it has been shown that "predictive power" is best for α < 0.1.H o wever, as detailed in [FRE 04], this approach suffers from known drawbacks. First, it depends on "magic values" for each weight often based on intuition or empirically determined. Second, it combines different scales of measurements that need to be normalized. Third, it combines non-commensurable criteria producing meaningless quantities. On the other hand, the lexicographic approach allows us to assign different priorities to different objectives in a qualitative fashion. Notably, in our context logic models providing high predictive power are significantly more relevant than the sizes of such models. Thus, the lexicographic approach is very convenient to cope with the multi-objective nature of our optimization problem. Yet another popular approach is to look for Pareto optimal models. However, this method will lead to a large number of models providing either none or very low predictive power. For example, consider the Boolean logic model (V, φ) with φ = ∅,i . e . t h eempty model. Such am o d e li st r i v i a l l yc o n s i s t e n tw i t ha n yi n p u tP K N(V, E, σ) while it minimizes the objective function size,i . e . Θ size ((V, φ)) = 0.T h e r e f o r e ,(V, φ) is Pareto optimal although it does not provide any valuable information. Similarly, one can show that many other (non-empty)modelswillbeParetooptimalaswellalthoughtheyprovide very low predictive power. Hence, Pareto optimality is not well suited for our problem. Notwithstanding, other multi-objective optimization methods (cf. [MAR 04]) could be investigated in the future. To conclude, our lexicographic multi-objective optimization consists of minimizing first Θ rss ,andthenwithlowerpriorityΘ size :

(V, φ opt ) = arg min (V,φ)∈M (V,E,σ) (Θ rss ((V, φ), ξ), Θ size ((V, φ))).
(1.4)

Information provided by high-throughput data is intrinsically uncertain due to experimental errors. Therefore, one is not only interested in optimal models but in nearly optimal models as well. In this context, authors in [SAE 09] have considered Boolean logic models minimizing Θ rss within certain tolerance, e.g. 10% of the minimum residual. Next, they argue that all models found can explain the data similarly or equally well if one take into account the experimental error. Notice that, in the aforecited work the optimization is addressed using a genetic algorithm. Hence, "minimum residual" refers to the minimum found during the execution of the algorithm which is not necessarily the global minimum. Moreover, due to the incompleteness of stochastic search methods, it is very likely that certain solutions within the allowed tolerance are not found. In practice, one can execute the genetic algorithm several times in order to overcome this issue to some extent. Nonetheless, as we show later in this chapter, a significant number of models may bemissingevenafterseveralexecutions. Similarly but in the context of quantitative modeling (based on ordinary differential equations) and using a simulated annealing algorithm, authors in [CHE 09] have elaborated upon the same argument. Interestingly, despite the fact that the model appears to be non-identifiable in both contexts, viz. qualitative and quantitative modeling, biologically relevant insights have been reported in the two aforecited studies. Notably, minimization over size in (1.4) is based onOccam'srazorprinciple.Ontheonehand, one can consider that larger logic models overfit the available dataset by introducing excessive complexity [SAE 09, PRI 11]. On the other hand, one can argue that it is actually necessary to consider such "spurious" links in order to capture cellular robustness and complexity [STE 04]. Therefore, let (V, φ opt ) be a Boolean logical model as defined in (1.4). Then, considering that tolerance over residual and size may yield biologically relevant models, we are particularly interested in enumerating all (nearly) optimal Boolean logic models (V, φ) such that, Θ rss ((V, φ), ξ) ≤ Θ rss ((V, φ opt ), ξ)+t rss Θ size ((V, φ)) ≤ Θ size ((V, φ opt ))+t size with t rss and t size denoting the tolerance over residual and size, respectively.

Next, we introduce the notion of logical input-output behaviors. In practice, the enumeration of (nearly) optimal models often leads to a large number of logical networks, namely, (V, φ j ) with j =1,...,mand m 1.N o t a b l y ,e a c hφ j is a different mapping from variables to propositional formulas. However, it may happen (and it often happens) that for all v ∈ V R ,severallogicalnetworksdescribeexactlythesame response to every possible experimental condition (clamping assignments over variables V S ∪ V K ). In such a case, we say that those logical networks describe the same input-output behavior. To be more precise, recall that we consider a PKN (V, E, σ). Notice that in each experimental condition over (V, E, σ),everystimulusv ∈ V S and inhibitor v ∈ V K ,canbeeitherclampedornot.Thus,letusdenotewithC the space of all possible clamping assignments or experimental conditions C over (V, E, σ). Notably, the number of possible clamping assignments is given by |C| =2 |VS|+|VK | .

Then, let (V, φ j ), (V, φ j ′ ) be two (nearly) optimal Boolean logic models. Furthermore, let F j C and F j ′ C be the fixpoints of Ω (V,φj|C ) and Ω (V,φ j ′ |C ) reachable from A f , respectively. We say that (V, φ j ) and (V, φ j ′ ) describe the same logical input-output behavior if and only if F j C (v)=F j ′ C (v) for all v ∈ V R and C ∈ C.I m p o r t a n t l y ,t h i s abstraction allows us to group logical networks regardless of their "internal wirings" and focus on their input-output predictions. In practice, this also facilitates the analysis and interpretation of results whereas it provides a way to extract robust insights despite the high variability.

Encoding

In order to express and solve the multi-objective optimization described in (1.4) by using ASP, one needs to discretize the function defined in (1.2). A very simple approach converts numerical data into binary data according to a threshold. Furthermore, we propose a finer multi-valued discretizationscheme. Inf act,theonlynon-inte ger variables in (1.2) are the experimental observations P Ci (v).T h e n ,w ea p p r o x i m a t e these values up to 1 10 k introducing a parametrized approximation function δ k (e.g. using the floor or closest integer functions). Next, we define the discrete residual Θ rss k as

Θ rss k ((V, φ), ξ)= n i=1 v∈dom(PC i ) 10 k δ k (P Ci (v)) -10 k π i (v) 2 .
(1.5)

The minimizations of Θ rss and Θ rss k may yield different Boolean logic models. Nonetheless, one can prove that finding all models minimizing Θ rss k within a certain tolerance allows us to find all models minimizing Θ rss as well.

Let (V, E, σ) be a PKN and let (V, H) be the directed hypergraph expanded from it. Recall that with P(v) we denote the powerset of the signed predecessors of v ∈ V , namely, 2 Pred(v) .W e r e p r e s e n t t h e d i r e c t e d h y p e r g r a p h (V, H) using predicates node/2, hyper/3,andedge/3.T h ef a c t snode(v,s P(v) ) map nodes v ∈ V to their corresponding sets of signed predecessors P(v),f a c t shyper(s P(v) ,s p ,l) associate P(v) with its sets p ∈ P(v) where l denotes their cardinalities, facts edge(s p ,v,1) associate the set p with (v, 1) ∈ p,a n df a c t sedge(s p ,v,-1) associate the set p with (v, -1) ∈ p.N o t et h a te a c hs (•) stands for some arbitrary but unique name in its respective context here. Facts over predicates stimulus/1, inhibitor/1,a n d readout/1 denote nodes in V S , V K ,a n dV R respectively. Next, let ξ =( C i ,P Ci ) be an experimental dataset over (V, E, σ) with i =1 ,...,n.R e c a l l t h a t e a c h C i is a clamping assignment over variables in V S ∪ V K .T h e n , w e e x t e n d o u r r e presentation of clamping assignments given before in order to consider several experimental conditions simultaneously. Towards this end, we represent experimental conditions as facts over predicate clamped/3,namelyclamped(i,v,C i (v)) for all v ∈ dom(C i ) and i =1 ,...,n. Finally, let k define the discretization scheme. We represent discretized experimental observations as facts over predicate obs/3,namely, obs(i,v,10 k δ k (P Ci (v))) for all v ∈ dom(P Ci ) and i =1,...,n.W eu s et h ep r e d icate dfactor/1 to denote the discretization factor 10 k .

Using the discretization scheme provided by k =1 ,L i s t i n g1 . 7s h o w st h ei nstance representation for our toy example. That is, the (signed) directed hypergraph in Figure 1.2(b) and the dataset given in (1.1).

Listing 1.7: Toy example input instance (toy.lp)

node(e,1). node(d,2). node(g,3). node(f,4). node(a,5). node(b,6). node(c,7).

hyper(1,1,1). hyper(2,1,1). hyper(1,8,2). hyper(2,13,2). hyper(1,2,1). hyper(2,4,1). hyper(1,9,2). hyper(2,11,2). hyper(1,3,1). hyper(2,5,1). hyper(1,10,2). hyper(2,12,2). hyper(3,5,1). hyper(4,6,1). hyper(3,14,2). hyper(1,16,3). hyper(3,6,1). hyper(4,7,1). hyper(4,15,2). hyper(2,17,3). edge(1,b,1). edge(2,c,1). edge(3,g,-1). edge(4,a,1). edge(5,c,-1). edge(6,e,1). edge(7,d,1). edge(8,b,1). edge(8,c,1). edge(9,b,1). edge(9,g,-1). edge(10,c,1). edge(10,g,-1). edge(11,a,1). edge (11,b,1). edge(12,a,1). edge(12,c,-1). edge(13,b,1). edge(13,c,-1). edge(14,e,1). edge(14,c,-1). edge(15,d,1). edge(15,e,1). edge(16,b,1). edge(16,c,1). edge(16,g,-1). edge(17,a,1). edge(17,b,1). edge(17,c,-1).

clamped (1,a,1). clamped (1,c,1). clamped (2,a,1). clamped (2,c,1). clamped ( 2,d,). clamped (3,a,1). clamped (4,a,1). clamped (4,b,1). obs(1,f,9). obs(1,g,0). obs(2,f,1). obs(2,g,9). obs(3,f,0). obs(3,g,1). obs(4,f,10). obs(4,g,8).

stimulus (a). stimulus (b). stimulus (c). inhibitor (d). readout(f). readout (g). dfactor (10).

Next we describe our encoding for solving the learning of Boolean logic models as described in the previous section. Our ASP encoding is shown in Listing 1.8. Listing 1.8: Logic program for learning Boolean logic models (learning.lp) variable (V) :-node(V,_). formula (V,I) :-node(V,I); hyper(I,_,_). {dnf(I,J) : hyper(I,J,N)} :-formula(V,I). clause(J,V,S) :-edge(J,V,S); dnf(_,J).

path(U,V) :-formula (V,I); dnf(I,J); edge(J,U,_). path(U,V) :-path(U,W); path(W,V). :-path(V,V).

:-dnf(I,J); edge(J,V,_); not stimulus (V); not path(U,V) : stimulus (U). :-path(_,V); not readout (V); not path(V,U) : readout (U).

exp(E)

:-clamped (E,_,_). clamped (E,V,-1) :-exp(E); stimulus (V); not clamped (E,V,1). clamped (E,V)

:-clamped(E,V,_). free(E,V,I)

:-formula(V,I); dnf(I,_); exp(E); not clamped (E,V).

eval(E,V, S) :-clamped (E,V,S). eval(E,V, 1) :-free(E,V,I); eval(E,W,T) : edge(J,W,T); dnf(I,J). eval(E,V,-1) :-not eval(E,V,1); exp(E); variable (V).

rss(D,V, 1, (F-D)**2) :-obs(E,V,D); dfactor (F). rss(D,V,-1, D**2) :-obs(E,V,D).

#minimize{L@1, dnf,I,J : dnf(I,J), hyper(I,J,L)}. #minimize{W@2, rss,E,V : obs(E,V,D), eval(E,V,S), rss(D,V,S,W)}.

:-formula(V,I); hyper(I,J1,N); hyper(I,J2,M); N < M, dnf(I,J1); dnf(I,J2); edge(J2,U,S) : edge(J1,U,S).

:-formula(V,I); dnf(I,J); edge(J,U,S); edge(J,U,-S).

#const maxsize = -1. #const maxrss = -1.

:-maxsize >= 0; maxsize + 1 #sum {L,dnf,I,J : dnf(I,J), hyper(I,J,L)}.

:-maxrss >= 0; maxrss + 1 #sum {W,rss,E,V : obs(E,V,D), eval(E,V,S), rss(D,V,S,W)}.

#show formula/2. #show dnf/2. #show clause/3.

Lines 1-4 define rules generating the representation of a logical network as described in Section 1.3. Line 1 simply projects node names to the predicate variable/1. In Line 2 every node v ∈ V having non-zero indegree is mapped to a formula φ(v). Next, in Line 3 each set of signed predecessors p ∈ P(v) is interpreted as an abducible conjuctive clause in φ(v).T h e n ,i nL i n e4i fp ∈ P(v) has been abduced, predicates clause/3 are derived for every signed predecessor in p. 6 Let us illustrate this on our toy example. In order to describe the mappings e → b ∨ c and g → e ∧ ¬c, one would generate a candidate answer set with atoms dnf(1,1), dnf(1,2) and 6. Notice that predicates clause/3 are only used for the sake of interpretation. One could simply replace Line 45 at the end of the encoding with #show clause(J,V,S) : edge(J,V,S); dnf(_,J). and remove Line 4. dnf(3,14) (from Line 2 we derive formula(e,1) and formula(g,3)). Note that this also force to have atoms clause(1,b,1), clause(2,c,1), clause(14,e,1) and clause(14,c,-1).L i n e s6 -8e l i m i n a t ec a n d i d a t ea n s w e rs e t sd e s c r i b i n gl o g i c models with feedback-loops. Paths from u to v are represented over predicate path/2 and derived recursively. Thus, the integrity constraint in Line 8 avoids self-reachability in the Boolean logic models. Next, the constraint in Lines 9-10 ensures that for any variable u defined in φ all variables w ∈ φ(u) are reachable from some stimuli variable. Whereas the constraint in Line 11 guarantees that every variable u defined in φ reaches some readout variable. Notice that att hispoint,w eha v earepres entationof the search space of Boolean logic models M (V,E,σ) .

Lines 13-21 elaborate on the rules from Listings 1.3, 1.4 and 1.5 given in Section 1.3 in order to consider several clamping assignments simultaneously and compute the fixpoint for each of them accordingly. To be more precise, the response under each experimental condition isr e p r e s e n t e do v e rp r e d i c a t e seval/3,n a m e l y eval(i,v,s) for experimental condition C i if variable v is assigned to s.I nL i n e s 23-24 we compute the possible differences (square of residuals) between Boolean predictions and the corresponding experimental observations. We denote such differences over predicate rss/4,n a m e l yrss(o,v,t,r) for a residual r with respect to the experimental observation o if the Boolean prediction for v ∈ V is the truth value t ∈ {1, -1}.N o t et h a ts u c hp r e d i c a t e sa r ei n d e p e n d e n tf r o me v e r yc a n d i d a t ea n s w e r set, that is, they can be deduced during grouding. For our example, due to the experimental condition C 2 we have rss(1,f,1,81), rss(1,f,-1,1), rss(9,g,1,1) and rss(9,g,-1,81).Therefore,ifthefixpointforf under the experimental condition C 2 is 1,t h er e s i d u a lw o u l db e81,w h e r e a si ft h efi x p o i n ti s-1,t h er e s i d u a li s only 1.A n a l o g o u s l y ,b u ti nt h eo p p o s i t ew a yt h es a m eh o l d sf o rg.N e x t ,w ed e s c r i b e our lexicographic multi-objective optimization. In Line 26 we declare with lower priority (@1)theminimizationoverthesizeoflogicmodels(Eq.(1.3)). Meanwhile,in Lines 27 we declare, with higher priority (@2), the minimization of the residual sum of squares between the Boolean predictions and experimental observations (Eq. (1.5)).

Lines 29-32 define two relatively simple symmetry-breaking constraints which are particularly relevant during the enumeration of (nearly) optimal solutions. Essentially, these integrity constraints eliminate answer sets describing "trivially" equivalent Boolean logic models with respect to their logical input-output behavior. The constraint in Lines 29-30 eliminates solutions by checking inclusion between conjunctions. For example, for two variables v and w,theformulav ∨ (v ∧ w) is logically equivalent to v and hence, the latter is preferred. Next, the constraint in Line 31 simply avoids solutions having mappings in the Boolean logic models of the form v ∧ ¬v. Notably, other logical redundancies could be considered as well. However, a complete treatment of redundancies would lead to the NP-complete problem known as minimization of Boolean functions [MCC 56] Lines 34-41 define a rather "standard" mechanism in order to enumerate solutions within given boundaries. Lines 34-35 simply define two constants describing the boundaries for each optimization criterion which are by default set to -1.Lines37-38 define an integrity constraint in order to eliminate solutions describing Boolean logic models (V, φ) if maxsize ≥ 0 and maxsize+1 ≤ Θ size ((V, φ)).A n a l o g o u s l y ,L i n e s 40-41 define an integrity constraint in order to eliminate solutions describing Boolean logic models (V, φ) if maxrss ≥ 0 and maxrss +1≤ Θ rss ((V, φ), ξ).

Solving

In Listing 1.9 we show the optimum answer set found for the toy instance described in Listing 1.7. 7 In this case, the optimum answer set is the thirteenth answer set inspected by the solver (Answer: 13). Such answer set describes the Boolean logic model given in Fig. 1.2(c). Furthermore, the values for the optimization criteria are given ordered by their priorities (Optimization: 88 7). That is, 88 for the discretized residual sum of squares (Eq. (1.5)), and 7 for the model size (Eq. (1.3)).

Listing 1.9: Learning an optimum Boolean logic model $g r i n g ot o y . l pl e a r n i n g . l p|c l a s p--q u i e t = 1 clasp version 3.0.2 Reading from stdin Solving ... Answer: 13 formula (e,1) formula (d,2) formula(g,3) formula (f,4)\ dnf (1,1) dnf (1,2) dnf (2,4) dnf (3,14) dnf (4,15)\ clause(1,b,1) clause(2,c,1) clause(4,a,1)\ clause (14,e,1) clause (14,c,-1) clause (15,d,1) clause (15,e,1) Optimization: 88 7 OPTIMUM FOUND

Models

: 13 Optimum : yes Optimization : 88 7 Calls : 1 Time : 0.005s (Solving :0.00s 1st Model :0.00s Unsat :0.00s) CPU Time : 0.000s

Next, the enumeration capabilities of an ASP solver like clasp [GEB 07] can be used to find not only one optimal model but all (nearly) optimal models as described earlier. Considering tolerance t rss =8(∼ 10% of the optimum residual sum of squares) and size tolerance t size =3,weenumerateallmodelssuchthat

Θ rss k ((V, φ), ξ) ≤ Θ rss k ((V, φ opt ), ξ)+t rss =96 Θ size ((V, φ)) ≤ Θ size ((V, φ opt )) + t size =10.
7. Using the option --quiet=1 only the last (optimum) answer set is printed. Notice that the solver prints all the atoms in the answer set in a single line but we have (manually) introduced breaklines to improve readability.

In this example, there are 5 (nearly) optimal Boolean logic models as we show in Listing 1.10. 8 Interestingly, even for this small example, the symmetry-breaking constraints make a significant difference. We note that running the same program but without the symmetry-breaking constraints, yields 17 Boolean logic models instead of only 5. Notably, in real-world problem instances, exploiting these symmetries significantly reduces the number of solutions (without missing any input-output behavior) and hence, it facilitates their post processing and interpretation. Once we have enumerated all (nearly) optimal Boolean logic models with respect to certain tolerances, we can identify the logical input-output behaviors they describe. Towards this end, we have developed a simple algorithm that (usingASP)systematicallycomparesallpairs of models looking for at least one experimental condition, i.e. a clamping assignment, generating a different response over the readouts nodes. We refrain from showing here the algorithm and additional encoding. Nonetheless, it is worth noting that the ASP encoding for deciding whether two Boolean logic models, there exists at least one experimental condition generating a different response over the readout nodes, is a rather straightforward extension from the rules given in Section 1.3. Following with our example, over the 5 (nearly) optimal Boolean logic models enumerated in Listing 1.10, we found 3 logical input-output behaviors.

Listing 1.10: Enumeration of all (nearly) optimal Boolean logic models $g r i n g ot o y . l pl e a r n i n g . l p-cm a x r s s = 9 6-cm a x s i z e = 1 0| \ clasp --opt -mode=ignore -n0 --quiet clasp version 3.0.2 Reading from stdin Solving ... SATISFIABLE

Models

: 5 Calls : 1 Time : 0.002s (Solving : 0.00s 1st Model: 0.00s Unsat: 0.00s) CPU Time : 0.000s

Importantly, results on a real-case study [GUZ 13] underscore the importance of exploring exhaustively the family of models and take into account experimental noise in order to obtain an adequate picture of the feasible model solutions. Briefly, in the aforecited work it is shown that if the experimental error is considered, several thousands of Boolean logic models fit the available data similarly well. Nonetheless, such alargenumberofmodelscanbegroupedintolessthanahundredlogicalinput-output behaviors. Next, these behaviors have been characterized in terms of the number of Boolean logic models they gather and their fitness to data. Moreover, it was found that for 30% of the space of possible inputs, all behaviors agree on the given outputs.

8. Option --opt-mode=ignore tells the solver to ignore optimize statements; option -n0 tells the solver to enumerate all solutions; and option --quiet avoids printing enumerated solutions.

Hence, in practice this approach may provide a way to extract robust insights despite the high variability. Also, thanks to our exhaustive characterization of these models, we can determine unambiguously which logical interactions are functional in all or none of the models and determine groups of mutually-exclusive mechanisms. Hence, our approach permits the study of the internal combinatorics leading to the variability of the system functioning and provides a tool to suggest new experiments towards the discrimination of plausible input-output behaviors given the available experimental setup. On the computational side, from multiple independent runs (1000 runs with an average of 1000 seconds per run) of the existing genetic algorithm implementation [TER 12], only 20% of them have converged to Boolean logic models within the allowed tolerance. Furthermore, among these runs, the genetic algorithm has retrieved approximately half of the plausible input-output behaviors identified using ASP with an evident bias towards the most common ones. Hence, those behaviors described only by a few logical networks are very unlikely to be found with such stochastic approaches.

Minimal intervention strategies

Background Once a biological system has been properly characterized in agreement with prior knowledge and experimental observations, a major challenge in systems biology is how to systematically control its state, leading to the field of experimental design. Importantly, progress in this area may have a crucial impact on bio-medical research, drug target identification and diagnosis. In fact, the problem of identifying "keyplayers" in biological systems has been addressed for metabolic, gene, and signaling networks. However, the underlying mathematical formalisms for each of these biological networks allow for different computational approaches.

In the context of logic-based models for signaling networks, this question has been recently addressed by [ABD 08, SAM 10, WAN 11,LAY 11]. Among them, in what follows we focus on the problem defined in [SAM 10]. Based on earlier work [KLA 06a] on metabolic networks, the notion of minimal intervention sets was introduced and dedicated algorithms were developed to compute them. Intuitively, an "intervention set" represents a set of knock-ins and knock-outs.E x a m p l e sf o rk n o c kins are mutations leading to constitutively activated species or a continuous stimulation with external signals whereas knock-outs may correspond to gene knock-outs or inhibition of a certain species by various experimental tools such as small-molecule drugs, antibodies, or RNAi. In fact, our notion of clamping assignments introduced in Section 1.3 was originally motivated as an abstraction (closer to standard logic terminology) of intervention sets. Then, we aim at identifying possible intervention sets leading to a logical steady state describing a specific biological outcome of interest. Furthermore, the aforecited work propose the usage of a three-valued logic which, though not mentioned, it correspond exactly to Kleene's logic [KLE 50]. It is also observed that under such a three-valued logic, a unique logical steady state follows for any (clamped) logical network and initial state. Hence, in contrast to the previous section where we have restricted ourselves to logical networks under Boolean logic and without feedback-loops, in this chapter, we consider logical networks under Kleene's logic and without any further restrictions. Therefore, the method presented in this section is not only relevant for Boolean logic models of immediate-early response, but for logical networks in general.

Unfortunately, the dedicated algorithms presented in [SAM 10] are computationally demanding due to the highly combinatorial mechanisms in logical signaling networks. Therefore, they are limited to compute small intervention sets and fail to scale over large-scale networks. In general, multiple interventions are necessary to cope with robustness and cellular complexity [STE 04]. Moreover, authors in [SAM 10] have considered looking for interventions in a single logical network. However, as we have shown above, if the inherent experimental noise is considered there are many logical networks compatible with a given dataset of experimental observations. Thus, identified interventions should fulfill the desired goals in every possible logical network. Concretely, the mentioned limitations make it hard to prove that the identified solutions are biologically robust to small perturbations of the system or its environment. Thus, in order to overcome such limitations, more elaborate and more powerful computational methods are needed towards large-scale systems and robust solutions.

Problem

Next, we provide a formal characterization of intervention strategies in logical signaling networks based on the notions introduced in Section 1.3. Given a logical network, the aim of an intervention strategy is to identify an intervention set that leads to a steady state satisfying a given goal under some side constraints. In fact, the concepts of an intervention (I), goal (G), and side constraints (C)canbecapturedas partial two-valued assignments. Moreover, both intervention sets and side constraints are considered clamping assignments as defined in Section 1.3 To be more precise, given a logical network (V, φ),anintervention scenario is a pair (G, C) of partial twovalue assignments over V where C is considered also as a clamping assignment, and an intervention set is a clamping assignment I over a set of intervention variables X ⊆ V .R e c a l lt h a tf o rt r u t ha s s i g n m e n t sA, B we defined the composition of assignments

A • B =(A \ B) ∪ B where B = {v → s | v → s ∈ B} and t = f , f = t, u = u.
Let (V, φ) be a logical network, let (G, C) be an intervention scenario, and X ⊆ V be a set of intervention variables. An intervention set I over X is an intervention strategy for (G, C) with respect to (V, φ),ifforsomej ≥ 0,wehavethat For illustration, let us consider the toy logical network in Figure 1.1(a) along with the intervention scenarios (G 1 ,C 1 )=( {o 1 → f ,o 2 → t}, {i 1 → t}) and (G 2 ,C 2 )=( {a → t}, ∅).I nt h i se x a m p l e ,t h efi r s ts c e n a r i or e q u i r e st h ei n h i b i t i o n of o 1 together with the activation of o 2 ,gi venthati 1 is stimulated. Furthermore, the second scenario requires the activation of a without any additional side constraints. Next, the intervention set {b → f ,e → f ,f → t} where b and e are inhibited and f is stimulated, satisfies both scenarios yielding the two steady states, respectively:

Ω j (V,φ|C•I ) (A u )=Ω j+1 (V,φ|C•I ) (A u ) G ⊆ Ω j (V,φ|C•I ) (A u ) with A u = {v → u | v ∈ V }.I nw o r d s ,Ω j (V,φ|C•I ) (A u )
⎧ ⎨ ⎩ i 1 → t i 2 → u a → t b → f c → f d → f e → f f → t g → t o 1 → f o 2 → t ⎫ ⎬ ⎭ ⎧ ⎨ ⎩ i 1 → u i 2 → u a → t b → f c → f d → f e → f f → t g → t o 1 → f o 2 → t ⎫ ⎬ ⎭
Authors in [SAM 10], were particularly interested in enumerating all minimal (bounded) intervention strategies with respect to a single logical network. However, as we have illustrated above and other authors have shown by considering real-world networks and data [SAE 09, CHE 09, GUZ 13], it often happens that the model is nonidentifiable. Therefore, as one can argue that several logical networks can describe a given biological system equally or similarly well, identified intervention strategies should fulfill all intervention scenarios in every possible logical network. Towards this end, herein we extend the problem settings in order to consider a family of logical networks, for instance, resulting from the enumeration of (nearly) optimal Boolean logic models described above. Now, let us define further intervention strategies relying on a finite family (V, φ i ) i∈N of logical networks, a finite family (G j ,C j ) j∈J of intervention scenarios and k some positive integer.

-Amulti-scenario intervention strategy for (G j ,C j ) j∈J wrt (V, φ i ) i∈N is an intervention strategy for each (G j ,C j ) wrt (V, φ i ) for each j ∈ J and i ∈ N .

-Abounded intervention strategy for (G j ,C j ) j∈J wrt (V, φ i ) i∈N and k is a multiscenario intervention strategy for (G j ,C j ) j∈J wrt (V, φ i ) i∈N of cardinality k ′ ≤ k.

-Aminimal bounded intervention strategy for (G j ,C j ) j∈J wrt (V, φ i ) i∈N and k is a ⊆-minimal multi-scenario intervention strategy for (G j ,C j ) j∈J wrt (

V, φ i ) i∈N of cardinality k ′ ≤ k.
In what follows, we focus on the enumeration of all minimal (bounded) intervention strategies for given families of intervention scenarios (G j ,C j ) j∈J and logical networks (V, φ i ) i∈N .

Encoding

Again, the representation of the problem instance is an extension from the one described in Listing 1.1 in order to describe a finite family of logical networks and clamping assignments. To be more precise, instead of having facts over predicates formula/2,weconsiderfactsoverpredicatesformula/3 as follows. Let (V, φ i ) i∈N be a finite family of logical networks. The facts formula(i,v,s φi(v) ) map variables v ∈ V to their corresponding formulas φ i (v) for each i ∈ N .M e a n w h i l e ,f a c t so v e r predicates variable/1, dnf/2 and clause/3 remain the same as in Listing 1.1. We use facts over predicate candidate/1 to denote the intervention variables that can be part of an intervention set. This allows us to control on which variables interventions are permitted, for example one can exclude interventions over constrained or goal variables. Next, we represent the family of intervention scenarios (G j ,C j ) j∈J using predicates scenario/1, goal/3,andconstrained/3.Thefactsscenario(j) denote the scenarios to consider. The facts goal(j,v,s) with s =1(resp. s

= -1)if G j (v)=t (resp. G j (v)=f )andconstrained(j,v,s) with s =1(resp. s = -1) if C j (v)=t (resp. C j (v)=f

variable (i1). variable (i2). variable (o2). variable (o1). variable (a). variable (b). variable (c). variable (d). variable (e). variable (f). variable (g). candidate (i2). candidate (b). candidate (c). candidate (d). candidate (e). candidate (f). candidate (g).

formula (1,a,0). formula (1,b,2). formula (1,c,1). formula ( 1,d,4). formula (1,e,3). formula (1,f,6). formula (1,g,5). formula (1,o1,4). formula (1,o2,7). dnf (0,5). dnf(1,6). dnf (1,0). dnf (2,3). dnf(3,7). dnf (4,1). dnf(5,2). dnf (6,4). dnf (6,6). dnf(7,4). clause(0,b,1). clause (1,c,1). clause(2,f,1). clause(3,a,1). clause(3,i1,1). clause (4,g,1). clause(5,d,-1). clause(6,e,1). clause(7,i2,1). clause (7,i1,-1). scenario (1). goal(1,o1,-1). goal(1,o2 ,1). constrained (1,i1,1). scenario (2). goal(2,a,1).

Next we describe our encoding for solving the minimal intervention set problem as described earlier. Our ASP encoding is shown in Listing 1.12.

Listing 1.12: Logic program for finding intervention strategies (control.lp) goal(T,S)

:-goal(_,T,S). goal(T)

:-goal(T,_). constrained (Z,E) :-constrained (Z,E,_). constrained (E) :-constrained (_,E). model(M) :-formula(M,_,_). formula (W,D) :-formula(_,W,D).

satisfy (V,W,S) :-formula (W,D); dnf(D,C); clause(C,V,S). closure (V,T) :-goal(V,T). closure (V,S*T) :-closure (W,T); satisfy(V,W,S); not goal(V,-S*T). {i n t e r v e n t i o n ( V , S ):closure (V,S) , candidate (V) }. :-intervention(V,1); intervention(V,-1). intervention(V) :-intervention(V,S). eval(M,Z,V,S) :-scenario (Z); intervention(V,S); model(M). eval(M,Z,E,S) :-model(M); constrained (Z,E,S); not intervention(E). free(M,Z,V,D) :-formula(M,V,D); scenario (Z); not constrained (Z,V); not intervention(V). eval_clause (M,Z,C,-1) :-clause(C,V,S); eval(M,Z,V,-S); model(M). eval(M,Z,V, 1) :-free(M,Z,V,D); dnf(D,C); eval(M,Z,W,T) : clause(C,W,T). eval(M,Z,V,-1) :-free(M,Z,V,D); eval_clause (M,Z,C,-1) : dnf(D,C).

:-goal(Z,T,S); model(M); not eval(M,Z,T,S). #const maxsize=0. :-maxsize >0; maxsize + 1 { intervention(X) }.

#show intervention/2.

In Lines 1-6 we define auxiliary domain predicates used in the remainder of the encoding. Lines 8-10 deserve closer attention since they allow us to reduce significantly the search space of candidate solutions. We incorporate a preprocessing step introduced in [SAM 10] that prunes variable assignments that can never be part of a minimal intervention strategy. The idea is to inductively collect all assignments that could be used to support a goal. First we gather all assignments that make a literal in a clause true and associate it with variable of the associated DNF (Line 8). Starting from the assignments that can satisfy a goal literal directly (Line 9), we inductively consider variable assignments (Line 10) that can support the assignments collected so far. Let us illustrate this on our toy example. In order to satisfy goal(1,o2,1),onew oul dne v er consider to intervene variables f or g negatively. Since both reach o2 positively, only positive interventions on them could help. The same happens for variable e.Ho we ver , since e also reaches o1 positively and we have goal(1,o1,-1), a negative intervention of e could help for this goal. Next, we use a choice rule in Line 12 to generate candidate solutions. We only choose interventions collected in the preprocessing step above. The integrity constraint in Line 13 eliminates contradictory interventions, e.g. intervention(e,1) and intervention(e,-1).WhereasLine14simplyprojects the intervention set to the intervened variables regardless of their signature. For example, one could generate the intervention set consisting of intervention(e,1) and intervention(c,-1).

In lines 16-26 elaborate on the rules from Listings 1.3, 1.4 and 1.6 given in Section 1.3 in order to consider several logical networks simultaneously and compute the fixpoint for each of them accordingly. To be more precise, we need to describe which variables are clamped (in all networks) according to the side constraints C j in each scenario j and the intervention set I,namely ,(V, φ i | Cj •I ).T o w a r d st h i se n d ,w eu s e the predicate eval/4,namelyeval(i,j,v,s) to represent that in the network (V, φ i ) and intervention scenario (G j ,C j ) the variable v is clamped to value s.F ol l o w i ngt he previous example, this will generate predicates eval(1,1,i1,1), eval(1,1,e,1), eval(1,2,e,1), eval(1,1,c,-1) and eval(1,2,c,-1).Theremainingrulesare adapted accordingly in order to compute for each logical network (V, φ i ) and intervention scenario (G j ,C j ),t h ec o r r e s p o n d i n gfi x p o i n to fΩ (V,φi|C j •I ) .F o r o u r e xample, we can see how the positive intervention over e is propagated as follows. Since variable f is not intervened, its formula (f → e ∨ g)d e s c r i b e db yp r e d i c a t e s formula(1,f,6), dnf(6,4), dnf(6,6), clause(4,g,1) and clause(6,e,1),is "free" in all scenarios, namely, free(1,1,f,6) and free(1,2,f,6).Furthermore, given that we have dnf(6,6) related only to clause(6,e,1) and e was intervened positively, the truth value for f is propagated in all scenarios regardless of g,namely , eval(1,1,f,1) and eval(1,2,f,1).

Line 28 declares an integrity constraint in order to eliminate answer sets describing intervention sets that do not satisfy all goals in every logical network and intervention scenario. Finally, the statements in Line 30 and 31 allows us to optionally bound the problem by considering only intervention sets up to a given size.

Solving

Normally ASP solvers allow for computing cardinality-minimal solutions whereas we are interested in finding subset-minimal solutions. In [KAM 13] we have shown how one can overcome this limitation by means of meta-programming and disjunctive logic programs [GEB 11a] or by using a specialized solver like hclasp [GEB 13]. However, herein we leverage thef u n c t i o n a l i t yr e c e n t l yi n t r o d u c e di nclasp 3s e r i e s which allow for computing subset-minimal solutions out-of-the-box by incorporating the features from hclasp.Importantly ,thisrequirestouseveryspecificcommand-line options for clasp.W er e f e rt h er e a d e rt oclasp's documentation for more details.

In Listing 1.13 we show the intervention strategies found for the toy instance described in Listing 1.11.

Listing 1.13: Computing all MISs for the toy instance $g r i n g oc o n t r o l . l pt o y . l p| \ clasp --dom -pref=32 --dom -mod=6 --heu=domain -n0 clasp version 3.0.2 Reading from stdin Solving ... Answer: 1 intervention(e,-1) intervention(f,1) intervention(b,-1) Answer: 2 intervention(i2,-1) intervention(f,1) intervention(b,-1) Answer: 3 intervention(e,-1) intervention(g,1) intervention(b,-1) Answer: 4 intervention(i2,-1) intervention(g,1) intervention(b,-1) Answer: 5 intervention(d,-1) intervention(f,1) intervention(b,-1) Answer: 6 intervention(g,1) intervention(d,-1) intervention(b,-1) Answer: 7 intervention(c,-1) intervention(f,1) Answer: 8 intervention(g,1) intervention(c,-1) Answer: 9 intervention(e,1) intervention(c,-1) SATISFIABLE

Models

: 9 Calls : 1 Time : 0.001s (Solving : 0.00s 1st Model: 0.00s Unsat: 0.00s) CPU Time : 0.000s

Interestingly, experiments over real-world biological networks in [KAM 13] have shown that our approach outperforms the previous dedicated algorithms [SAM 10] in up to four orders of magnitude (for small number of interventions (≤ 3)stillfeasible for the algorithm). This was not very surprising since such algorithms are based on astandardbreadth-firstsearchusingadditionaltechniquesforsearchspacereduction. More importantly, using ASP we are able to search for significantly larger intervention strategies or even solve the unbounded problem (that is, no limit in the number of interventions). While considering a small number of interventions the number of solutions (that is, intervention strategies) is in the order of tens, with a larger number of interventions we have found thousands of feasible solutions. Furthermore, being able to solve the unbounded problem allows us to completely characterize the set of popular scripting language, e.g. python. 9 Analogously, the resulting answer sets from clasp can be converted back to "standard" formats for subsequent analysis with available tools or even visualized in order to facilitate their interpretations. Hence, towards this end we have implemented the python package caspo which is freely available for download. 10 The aim of caspo is to implement a pipeline for automated reasoning on logical signaling networks providing a powerful and easy-to-use software tool for systems biologists. In particular, both problems presented herein together with other related features are available. For the sake of illustration, in Listing 1.14 we show the usage of caspo over a real-world instance from the "DREAM" challenge [ PRI 11]. For more details on the installation, usage, and available features we refer the reader to the online documentation. More broadly, caspo is part of BioASP, a collection of python packages leveraging the computational power of ASP for systems biology. 11

Conclusion

Systems biology is a research field at the crossover of biology, informatics, and mathematics. Its central task is to decipher the functioning of the so-called biological networks mediating between DNA products, the environment, and the phenotype. Among these networks, signal transduction networks are crucial for the understanding of the response to extra-and intracellular perturbations. Notably, in order to cope with the increasing complexity of large-scale networks, the development of mathematical and computational modeling approaches is a major goal in the field. In this context, despite their high-level of abstraction, logical networks provide a powerful qualitative approach for modeling large-scale biological systems. Importantly, several authors have shown that the response of signaling networks can be appropriately modeled with (Boolean) logical networks. Nevertheless, several challenging problems remain open. For example, the question of identifying the precise logic rules underlying the system being studied or finding appropriate interventions allowing to control it. Existing approaches implementingdedicatedalgorithmstoaddresstheseproblemshave been relatively successful. However, they possess certain shortcomings regarding the scalability and exhaustiveness over large search spaces. In fact, this incompleteness in the search may significantly compromise the robustness of solutions found while it limits the insights provided to biologists.

In contrast to previous approaches, herein we propose to look for robust insights by reasoning over the complete space of feasible solutions. Towards this end, we strongly rely on methods from the area of knowledge representation and reasoning such as Answer Set Programming (ASP). Available systems nowadays provide a rich yet simple modeling language,h i g h -p e r f o r m a n c es o l v i n gc a p a c i t i e s ,a n da u t o m a t e d reasoning modes. Thus, ASP provides a powerful framework for addressing hard combinatorial problems in systems biology by reasoning over the complete search space. In order to illustrate this, we have shown how the aforementioned problems on logical networks can be successfully modeled and solved using ASP. Interestingly, in both cases the computational performance is significantly improved with respect to dedicated algorithms. But more importantly, the exhaustive nature of ASP allows us to find feasible solutions that were missing when using the existing methods. Altogether, our work constitutes a pipeline for automated reasoning on logical signaling networks providing robust insights to systems biologists.

Finally, as mentioned earlier, integrative modeling approaches considering multiple levels and time-scales of causation pose a very challenging goal in systems biology. In order to achieve this goal, we believe that more sophisticated computational methods are needed in order to integrate qualitative and quantitative knowledge. In particular, hybrid reasoning systems leveraging the expressiveness of several technologies and modeling approaches appears as a very promising track for future research and development. Hopefully, advances on this subject will foster the usage of knowledge representation and reasoning methodologies in systems biology towards a better understanding of life.

Acknowledgements

We would like to thank our collaborators: Axel von Kamp, Carito Guziolowski, Federica Eduati, Jacques Nicolas, Julio Saez-Rodriguez, Martin Gebser, Regina Samaga, Roland Kaminski, Steffen Klamt, Sven Thiele, and Thomas Cokelaer. This work was partially funded by the projects ANR-10-BLANC-0218 and DFG grant SCHA 550/10-1.

Bibliography

[ABD 08] ABDI A., TAHOORI M. B., EMAMIAN E. S., "Fault diagnosis engineering of digital circuits can identify vulnerable molecules in complex cellular pathways", Science Signaling,vol.1,num.42,2008.

[AKU 00] AKUTSU T., MIYANO S., KUHARA S., "Inferring qualitative relations in genetic networks and metabolic pathways", Bioinformatics,vol.16,num.8,p.727-734,July 2000.

[ALE 10] ALEXOPOULOS L. G., SAEZ-RODRIGUEZ J., COSGROVE B., LAUFFENBURGER D. A., SORGER P., "Networks inferred from biochemical data reveal profound differences in toll-like receptor and inflammatory signaling between normal and transformed hepatocytes", Molecular & Cellular Proteomics,vol.9,num.9,p.1849-1865,2010.

[

1. 1 .

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2. Answer Set Programming at a glance . . . . . . . . . . . . . . . . . . . 1.3. Learn and control logical networks with ASP . . . . . . . . . . . . . . . 1.3.1. Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2. Reasoning on the response of logical networks . . . . . . . . . . . 1.3.3. Learning models of immediate-early response . . . . . . . . . . . 1.3.4. Minimal intervention strategies . . . . . . . . . . . . . . . . . . . . 1.3.5. Software toolbox: caspo . . . . . . . . . . . . . . . . . . . . . . . . 1.4. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5. Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6. Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Index .......................................... 5 Chapter 1

Figure 1 . 1 .

 11 Figure 1.1. Graphical representaion for logical networks and interaction graphs. Vertices represent biological entities and the interactions among them are represented as follows. Positive interactions are represented by an arrow (→)whereasnegativeinteractionsarerepresentedby aT -shape(⊣). (a) Exemplary logical network represented as a directed hypergraph. Directed hyperedges describe logical interactions. (b) Interaction graph underlying the logical network in (a). Directed edges describe causal interactions.

  and initial assignment Au or A f . corresponding evaluation. Finally, we use the predicate free/2 to represent the fact that every variable not clamped in C,i ss u b j e c tt ot h ec o r r e s p o n d i n gm a p p i n gφ(v). Both rules are shown in Listing 1.3. Listing 1.3: Clamped and free variables eval(V,S) :-clamped(V,S). free(V,D) :-formula(V,D); dnf(D,_); not clamped (V,_).

Figure 1 . 2 .

 12 Figure1.2. The green and red edges correspond to activations and inhibitions, respectively. Green nodes represent ligands that can be experimentally stimulated. Red nodes represent species that can be inhibited by using a drug. Blue nodes represent species that can be measured by using an antibody. White nodes are neither measured, nor manipulated. (a) At o y interaction (directed and signed) graph describing causal interactions. (b) Hypergraph expansion describing all plausible logical interactions based on the prior knowledge network in (a). (c) AB o o l e a nl o g i cm o d e ld e r i v e df r o mt h ep r i o rk n o w l e d g en e t w o r ki n( a )d e s c r i b i n g functional relationship defined by the mapping {d → a; e → b ∨ c; f → d ∧ e; g → e ∧ ¬c}.

  is a steady state of the clamped network (V, φ| C•I ) satisfying the goal conditions in G.N o t i c et h ec o m p o s ition of C • I indicating that clampings in the intervention set I overwrite clampings in the side constraints C.F i n a l l y ,t h eintervention set problem consists in deciding whether there is an intervention strategy for an intervention scenario (G, C) wrt a logical network (V, φ).

  )d e n o t et h er e s p e c t i v ei n t e r v e n t i o ng o a l sa n ds i d e constraints in each scenario (G j ,C j ).Listing 1.11 shows the instance representation of our toy example logical network in Figure1.1(a) together with the two intervention scenarios (G 1 ,C 1 )=( {o 1 → f ,o 2 → t}, {i 1 → t}) and (G 2 ,C 2 )=( {a → t}, ∅).N o t a b l y ,f o rt h es a k eo fu nderstanding, we consider a toy example with a single logical network. But in general, this instance representation and the logic program given below, support several logical networks.Listing 1.11: Toy example problem instance (toy.lp)

problem consisting of reconstructing all possible networks consis- tent with experimental time series data [DUR 11]. Gebser et al. have addressed the problem consisting of detecting inconsistencies and repairing in large biological net- works [GEB 11b, GEB 10]. Fayruzov et al. have used ASP to represent the dynamics in Boolean networks and find their attractors [FAY 11]. Ray et al. have integrated numerical and logical information in order to find the most likely states of

  Given a finite set V of propositional variables, we form propositional formulas from V with the connectives ⊥, ⊤, ¬, ∨,and∧ in the standard way. Further, we consider (partial) truth assignments over V mapping formulas to truth values {t, f , u} according to Kleene's semantics [KLE 50]. Clearly, two-valued assignments are restricted to range {t, f } according to classical (Boolean)l o g i cs e m a n t i c s . W er e c a l l the truth tables for classical (Boolean) and Kleene's logics in Table1.1 and Table1.2, respectively.

	a ∧ b a t t f b t f f f f	a ∨ b a t t t b t f f t f	a ¬a t f f t
	Table 1.1. Truth tables for classical (Boolean) logic.
	1.3. Learn and control logical networks with ASP
	1.3.1. Preliminaries
	Propositional logic and mathematical notation
	al. have proposed applying knowledge representation
	and reasoning methodologies to the problem of representing and reasoning about sig-
	naling networks [BAR 04]. More recently, several authors have addressed the ques-
	tion of pruning or identification of biological networks using ASP. Durzinsky et al.
	have studied the a biolog-
	ical system under various constraints [RAY 12]. Furthermore, Ray et al. have used
	an ASP system to propose revisions to metabolic networks [RAY 10]. Papatheodorou
	et al. have used ASP to integrate RNA expression with signaling pathway informa-
	tion and infer how mutations affect ageing [PAP 12]. Finally, Schaub and Thiele have
	first investigated the metabolic network expansion problem with ASP [SCH 09b] and
	recently, their work has been extendend and applied in a real-case study by Collet et
	al. [COL 13]. Altogether, this series of contributions illustrates the potential of ASP
	to address combinatorial search and optimization problems appearing in the field.
	Nonetheless, its strictly discrete nature poses interesting challenges for future work
	towards hybrid reasoning system allowing for qualitative and quantitative modeling.

Table 1 .

 1 2. Truth tables for Kleene's logic.

  APT 82] APT K. R., EMDEN M. H. V., "Contributions to the Theory of Logic Programming", ACM,vol.29,num.3,p.841-862,ACM,July 1982.

		Index
	answer set, 16	intervention goal, 40
	Answer Set Programming, 16	intervention set problem, 41
	ASP encoding, 23, 33, 42	intervention strategy, 40
	ASP solving, 37, 44 biological network, 12 Boolean model accuracy, 31 Boolean model accuracy (discrete), 33 Boolean model complexity, 31 Boolean network, 14, 20 Boolean prediction, 28 caspo, 46 clamping, 21, 24, 28, 36, 39, 40 clasp, 16, 37, 45 experimental condition, 28 experimental design, 39 feedback-loops, 22, 26, 29, 40 fixpoint, 22, 26, 29, 36, 44 genetic algorithm, 15, 32 gringo 4, 17 hclasp, 44 hypergraph, 20 immediate-early response, 23, 26, 40 inference, 15 input-output behavior, 32, 33, 38	key-players, 15, 39 knock-in, 27, 39 knock-outs, 27, 39 learning, 14, 27 logical network, 14, 19, 29, 40 multi-objective function for training Boolean models, 31 nearly optimal models, 32 phosphorylation assays, 13, 30 prior knowledge network, 26 readouts, 27 reasoning mode, 16, 48 reverse engineering, 14 search space, 27, 30, 43 semantics, 19, 40 signaling network, 13, 39 single-step operator for logic program, 22 steady state, 20, 22, 26, 39, 41 stimuli, 26, 39 systems biology, 11, 18, 47
	interaction graph, 21, 26	three-valued logic, 19, 22, 25, 40
	intervention, 15, 39	tolerance, 15, 32
	intervention constraints, 40	two-valued logic, 19, 22, 25

http://www.mat.unical.it/aspcomp2013/ASPStandardization

Also known as sum-of-products [KLA 06b].

More generally, a directed hyperedge is a pair (S, T ) with T ⊆ V .W ec o n s i d e rt h ep a r t i c u l a rc a s e where T is a singleton. Directed hypergraphs are sometimes referred to as "AND/OR graphs"[GAL 93]

The interested reader may notice the resemblance tosi n g l e -st e po p e r a t o r sf o rl o g i cp r o g r a msi n t r oduced in [APT 82] and [FIT 85] for two-and three-valued assignments respectively.

solutions. Notably, in the light of such a large number of solutions, the way to select among them arises. Herein, we have introduced an extension of the intervention sets problem in order to consider a family of plausible logical networks as those identified in the previous section. This way, we expect to reduce the number of solutions by selecting the more robust of them. That is, intervention sets satisfying each scenario in all logical networks. In practice, interactions graphs, experimental datasets, logical networks, and similar knowledge in systems biology is often (publicly) available in different kind of "standard" formats. Clearly, converting such a knowledge from their corresponding format to a set of logic facts, e.g. in Listing 1.7 and Listing 1.11, is a tedious and error-prone task if do it by hand. In fact, this can be easily automated by using any