Target identification, one of the steps of drug discovery, aims at identifying biomolecules whose function should be therapeutically altered in order to cure the considered pathology. This work proposes an algorithm for in silico target identification using Boolean network attractors. It assumes that attractors of dynamical systems, such as Boolean networks, correspond to phenotypes produced by the modeled biological system. Under this assumption, and given a Boolean network modeling a pathophysiology, the algorithm identifies target combinations able to remove attractors associated with pathological phenotypes. It is tested on a Boolean model of the mammalian cell cycle bearing a constitutive inactivation of the retinoblastoma protein, as seen in cancers, and its applications are illustrated on a Boolean model of Fanconi anemia. The results show that the algorithm returns target combinations able to remove attractors associated with pathological phenotypes and then succeeds in performing the proposed in silico target identification. However, as with any in silico evidence, there is a bridge to cross between theory and practice, thus requiring it to be used in combination with wet lab experiments. Nevertheless, it is expected that the algorithm is of interest for target identification, notably by exploiting the inexpensiveness and predictive power of computational approaches to optimize the efficiency of costly wet lab experiments.

Introduction

Drug discovery, as its name indicates, aims at discovering new drugs against diseases. This process can be segmented into three steps: i) disease model provision, where experimental models are developed, ii) target identification, where therapeutic targets are proposed, and iii) target validation, where the proposed therapeutic targets are assessed. This work focuses on the second step of drug discovery: target identification [START_REF] Knowles | Target selection in drug discovery[END_REF][START_REF] Mark | Target discovery[END_REF].

Given an organism suffering from a disease, target identification aims at finding where to act among its multitude of biomolecules in order to alleviate, or ultimately cure, the physiological consequences of the disease. These biomolecules on which perturbations should be applied are called targets and are targeted by drugs [START_REF] Imming | Drugs, their targets and the nature and number of drug targets[END_REF]. This raises two questions: which target should be therapeutically perturbed and what type of perturbation should be applied on it. Broadly, the functional perturbation of a target by a drug can be either activating or inactivating, regardless the way the drug achieves it.

One solution is to test all, or at least a large number of, biomolecules for activation and inactivation. Knowing that targeting several biomolecules is potentially more effective [START_REF] Anighoro | Polypharmacology: challenges and opportunities in drug discovery[END_REF][START_REF] Grant R Zimmermann | Multi-target therapeutics: when the whole is greater than the sum of the parts[END_REF], the number of possibilities is consequently huge. This rather brute-force screening can be refined with knowledge about the pathophysiology of interest by identifying potential targets based on the role they play in it [START_REF] Jackson | Mechanism-based target identification and drug discovery in cancer research[END_REF]. Even with this knowledge, experimentally assessing the selected potential targets through wet lab experiments is far from straightforward since such experiments are costly in time and resources [START_REF] Ki Kaitin | Deconstructing the drug development process: the new face of innovation[END_REF]. Fortunately, owing to their integrative power and low cost compared to wet lab experiments, in silico approaches appear as valuable tools in improving the efficiency of target identification [START_REF] Tang | Network pharmacology strategies toward multi-target anticancer therapies: from computational models to experimental design principles[END_REF][START_REF] Mak | Anti-cancer drug development: computational strategies to identify and target proteins involved in cancer metabolism[END_REF][START_REF] Sarker | In silico systems biology approaches for the identification of antimicrobial targets[END_REF][START_REF] Yao | In silico search for drug targets of natural compounds[END_REF][START_REF] Praveen | Targeting drug transporters-combining in silico and in vitro approaches to predict in vivo[END_REF][START_REF] Hua | In silico approaches to multi-target drug discovery[END_REF][START_REF] Vujasinovic | In silico dynamic molecular interaction networks for the discovery of new therapeutic targets[END_REF][START_REF] Chandra | Computational systems approach for drug target discovery[END_REF][START_REF] Saidani | Potential and limits of in silico target discovery-case study of the search for new antimalarial chemotherapeutic targets[END_REF][START_REF] Ulrik | Using computational modeling to drive the development of targeted therapeutics[END_REF][START_REF] Duckworth | In silico identification of novel therapeutic targets[END_REF][START_REF] Noble | Biological simulations in drug discovery[END_REF], as demonstrated through several works using various computational methods [START_REF] Carels | A computational strategy to select optimized protein targets for drug development toward the control of cancer diseases[END_REF][START_REF] Frangou | Molecular profiling and computational network analysis of taz-mediated mammary tumorigenesis identifies actionable therapeutic targets[END_REF][START_REF] Li | Computer-aided targeting of the pi3k/akt/mtor pathway: toxicity reduction and therapeutic opportunities[END_REF][START_REF] Nicklas | In silico identification of potential therapeutic targets in the tgf-β signal transduction pathway[END_REF][START_REF] Rao | Identification of novel drug targets in hpb38, hpp12, hpg27, hpshi470, hpsjm180 strains of helicobacter pylori: an in silico approach for therapeutic intervention[END_REF][START_REF] Bilachi S Ravindranath | In silico analyses of metabolic pathway and protein interaction network for identification of next gen therapeutic targets in chlamydophila pneumoniae[END_REF][START_REF] Wu-Lung R Yang | In silico drug screening and potential target identification for hepatocellular carcinoma using support vector machines based on drug screening result[END_REF][START_REF] Iadevaia | Identification of optimal drug combinations targeting cellular networks: integrating phospho-proteomics and computational network analysis[END_REF][START_REF] Betul | Drug target identification in sphingolipid metabolism by computational systems biology tools: metabolic control analysis and metabolic pathway analysis[END_REF][START_REF] On Koborova | In silico method for identification of promising anticancer drug targets[END_REF][START_REF] Madhukar S Dasika | A computational framework for the topological analysis and targeted disruption of signal transduction networks[END_REF].

However, the stumbling block of in silico approaches is that they are built from the available knowledge: not all is known about everything. Nevertheless, an impressive and ever increasing amount of biological knowledge is already available in the scientific literature, databases and knowledge bases such as, to name a few, DrugBank [START_REF] David S Wishart | Drugbank: a knowledgebase for drugs, drug actions and drug targets[END_REF], KEGG [START_REF] Kanehisa | Kegg: kyoto encyclopedia of genes and genomes[END_REF], PharmGKB [START_REF] Whirl-Carrillo | Pharmacogenomics knowledge for personalized medicine[END_REF], Reactome [START_REF] Croft | Reactome: a database of reactions, pathways and biological processes[END_REF] and TTD [START_REF] Chen | Ttd: therapeutic target database[END_REF]. In addition to the difficulty of integrating an increasing body of knowledge comes the inherent complexity of biological systems themselves [START_REF] Kitano | Systems biology: a brief overview[END_REF]: this is where computational tools can help owing to their integrative power [START_REF] Boissel | Modeling and medical product r& d[END_REF][START_REF] Boissel | Modelling methodology in physiopathology[END_REF][START_REF] Kitano | Computational systems biology[END_REF]. This interplay between wet lab and computational biology is synergistic rather than competitive [START_REF] Ventura | From in vivo to in silico biology and back[END_REF]. Since wet lab experiments produce factual results, they can be considered as trustworthy sources of knowledge. Once these factual pieces of knowledge are obtained, computational tools can help to integrate them and infer new ones. This computationally obtained knowledge can be subsequently used to direct further wet lab experiments, thus mutually potentiating the whole.

The goal of this work is to propose a computational methodology implemented in an algorithm for in silico therapeutic target discovery using Boolean network attractors. It assumes that Boolean network attractors correspond to phenotypes produced by the modeled biological network, an assumption successfully applied in several works [START_REF] Von Der Heyde | Boolean erbb network reconstructions and perturbation simulations reveal individual drug response in different breast cancer cell lines[END_REF][START_REF] Herman | Boolean network model for cancer pathways: predicting carcinogenesis and targeted therapy outcomes[END_REF][START_REF] Grieco | Integrative modelling of the influence of mapk network on cancer cell fate decision[END_REF][START_REF] Andres Moreno-Ramos | Modeling of the hypothalamic-pituitary-adrenal axis-mediated interaction between the serotonin regulation pathway and the stress response using a boolean approximation: a novel study of depression[END_REF][START_REF] Pau Creixell | Navigating cancer network attractors for tumor-specific therapy[END_REF][START_REF] Rodriguez | A boolean network model of the fa/brca pathway[END_REF][START_REF] Singh | Boolean approach to signalling pathway modelling in hgfinduced keratinocyte migration[END_REF][START_REF] Baverstock | A comparison of two cell regulatory models entailing high dimensional attractors representing phenotype[END_REF][START_REF] Naldi | Diversity and plasticity of th cell types predicted from regulatory network modelling[END_REF][START_REF] Thakar | Boolean models of within-host immune interactions[END_REF][START_REF] Ge | Boolean network approach to negative feedback loops of the p53 pathways: synchronized dynamics and stochastic limit cycles[END_REF][START_REF] Sahin | Modeling erbb receptor-regulated g1/s transition to find novel targets for de novo trastuzumab resistance[END_REF][START_REF] Maria | Boolean network model predicts cell cycle sequence of fission yeast[END_REF][START_REF] Kervizic | Dynamical modeling of the cholesterol regulatory pathway with boolean networks[END_REF][START_REF] Faure | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF][START_REF] Mendoza | A network model for the control of the differentiation process in th cells[END_REF][START_REF] Huang | Shape-dependent control of cell growth, differentiation, and apoptosis: switching between attractors in cell regulatory networks[END_REF]. Assuming that a phenotype is an observable state, and thus relatively stable, of a biological system and assuming that the state of a biological system results from its dynamics, a phenotype is likely to correspond to an attractor. This assumption can be stated for any dynamical model but, in this work, only Boolean networks are considered. Reasons are that, in their most basic form, Boolean networks do not require quantitative information [START_REF] Michelle L Wynn | Logic-based models in systems biology: a predictive and parameter-free network analysis method[END_REF] and that quantitative information is often not easy to obtain due to experimental limitations, particularly at the subcellular scale, the scale where drugs interact with their targets. Moreover, since synchronous Boolean networks are easier to compute than asynchronous ones [START_REF] Garg | Synchronous versus asynchronous modeling of gene regulatory networks[END_REF], this work only considers synchronous Boolean networks. This does not exclude the possibility, at a later stage, to extend the algorithm for both synchronous and asynchronous updating schemes.

For a biological network involved in a disease, two possible variants are considered: the physiological variant, exhibited by healthy organisms, which produces physiological phenotypes, and the pathological variant, exhibited by ill organisms, which produces pathological phenotypes or which fails to produce physiological ones. A physiological phenotype does not impair life quantity/ quality whereas a pathological phenotype does. It should be noted that the loss of a physiological phenotype is also a pathological condition. The physiological and pathological variants differ in that the latter results from the occurrence of some alterations known to be responsible for disorders. With a pathological variant, there are two non-exclusive pathological scenarios: pathological phenotypes are gained or physiological phenotypes are lost.

The primary goal of the proposed algorithm is to identify, in a pathological variant, target combinations together with the perturbations to apply on them, here called bullets, which render it unable to exhibit pathological phenotypes. The secondary goal is to classify the obtained bullets according to their ability at rendering the pathological variant able to exhibit the previously lost physiological phenotypes, if any.

It should be noted that this work fits into the encompassing field investigating how to control biological systems, a field with tremendous applications in biomedicine. Several endeavors based on qualitative modeling approaches have been made in this way [START_REF] Campbell | Stabilization of perturbed boolean network attractors through compensatory interactions[END_REF][START_REF] Qiu | On control of singleton attractors in multiple boolean networks: integer programming-based method[END_REF][START_REF] Srihari | Evolution and controllability of cancer networks: a boolean perspective[END_REF][START_REF] Chen | Finding optimal control policy in probabilistic boolean networks with hard constraints by using integer programming and dynamic programming[END_REF][START_REF] Kobayashi | Optimal control of gene regulatory networks with effectiveness of multiple drugs: A boolean network approach[END_REF][START_REF] Kobayashi | Symbolic approach to verification and control of deterministic/probabilistic boolean networks[END_REF], demonstrating its utility in investigating how to take control over pathologically disturbed biological systems.

Methods

This section introduces some basic principles, namely biological and Boolean networks, defines some concepts and then describes the proposed algorithm. An example network to illustrate how it works plus a case study to illustrate its intended applications are also described. Finally, details about implementation and code availability are mentioned.

Basic principles

Biological networks

A biological network is a way to conceptualize a set of interacting biological entities where entities are represented by nodes and interactions by edges [START_REF] Zhu | Getting connected: analysis and principles of biological networks[END_REF][START_REF] Barabasi | Network biology: understanding the cell's functional organization[END_REF]. It is based on graph theory [START_REF] Nikolaevich Bronshtein | Algorithms of graph theory[END_REF][START_REF] Huber | Graphs in molecular biology[END_REF][START_REF] Mason | Graph theory and networks in biology[END_REF], thus bringing formal tools to encode information about biological systems, particularly their topology [START_REF] Ma'ayan | Insights into the organization of biochemical regulatory networks using graph theory analyses[END_REF]. Moreover, being graphs, biological networks offer a convenient visualization [START_REF] Jill | Why a diagram is (sometimes) worth ten thousand words[END_REF] of the complex interconnections lying in biological systems. As said Napoleon Bonaparte: "A good sketch is better than a long speech." Mathematically, a network can be seen as a digraph G = (V, E) where V = {v 1 , . . . , v n } is the set of cardinality n containing exactly all the nodes v i of the network and where E = {(v i,1 , v j,1 ), . . . , (v i,m , v j,m )} ⊆ V 2 is the set of cardinality m containing exactly all the edges (v i , v j ) of the network. In practice, nodes represent entities and edges represent binary relations R ⊆ V 2 involving them: v i R v j . For example, in gene regulatory networks, nodes represent gene products and edges represent gene expression modulations [START_REF] Xiao | A tutorial on analysis and simulation of boolean gene regulatory network models[END_REF][START_REF] John | Gene expression networks[END_REF].

Boolean networks

While being conceptually simple, Boolean networks [START_REF] Saadatpour | Boolean modeling of biological regulatory networks: a methodology tutorial[END_REF] are able to predict and reproduce features of biological systems and then to bring relevant insights [START_REF] Albert | Boolean modeling: a logic-based dynamic approach for understanding signaling and regulatory networks and for making useful predictions[END_REF][START_REF] Wang | Boolean modeling in systems biology: an overview of methodology and applications[END_REF][START_REF] Albert | Boolean network simulations for life scientists[END_REF][START_REF] Bornholdt | Boolean network models of cellular regulation: prospects and limitations[END_REF][START_REF] Huang | Genomics, complexity and drug discovery: insights from boolean network models of cellular regulation[END_REF][START_REF] Leclerco | Boolean analysis of cell regulation networks[END_REF]. This makes them an attractive and efficient approach, especially when the complexity of biological systems renders quantitative approaches unfeasible due to the amount of quantitative details they require. As their name indicates, Boolean networks are based on Boolean logic [START_REF] Boole | The mathematical analysis of logic, being an essay towards a calculus of deductive logic[END_REF] and, like biological networks, are also based on graph theory: nodes represent Boolean variables and edges represent interdependencies between them.

Mathematically, a Boolean network is a network where nodes are Boolean variables x i and where edges (x i , x j ) represent the binary is input of relation:

x i is input of x j . Each x i has b i ∈ [[0, n]] inputs x i,1 , . . . , x i,bi
. The variables which are not inputs of x i have no direct influence on it. If b i = 0 then x i is a parameter and does not depend on other variables. At each iteration k ∈ [[k 0 , k end ]] of the simulation, the value x i (k) ∈ {0, 1} of each x i is updated to the value x i (k + 1) using a Boolean function f i and the values x i,1 (k), . . . , x i,bi (k) of its inputs, as in the following pseudocode:

1 for k ∈ [[k 0 , k end -1]] do 2 x 1 (k + 1) = f 1 (x 1,1 (k), . . . , x 1,b1 (k)) 3 . . . 4 x n (k + 1) = f n (x n,1 (k), . . . , x n,bn (k))
5 end for which can be written in a more concise form:

1 for k ∈ [[k 0 , k end -1]] do 2 x(k + 1) = f (x(k))
3 end for where f = (f 1 , . . . , f n ) is the Boolean transition function and x = (x 1 , . . . , x n ) is the state vector. The value x(k) = (x 1 (k), . . . , x n (k)) ∈ {0, 1} n of x at k belongs to the state space S = {0, 1} n which is the set of cardinality 2 n containing exactly all the possible states.

If the values of all the x i are updated simultaneously at each k then the network is synchronous, otherwise it is asynchronous. With synchronous Boolean networks, x(k) has a unique possible successor x(k + 1): synchronous Boolean networks are deterministic. In the particular case where k = k 0 , x(k 0 ) = x 0 is the initial state and, in deterministic dynamical systems, determines entirely the trajectory w = (x(k 0 ), . . . , x(k end )). In this work, it is assumed that k 0 = 1, so w is a sequence of length k end resulting from the iterative computation of x(k) from k 0 up to k end . This iterative computation can be seen as the discretization of a time interval: Boolean networks are discrete dynamical systems as they simulate discretely the time course of the state vector.

The set A = {a 1 , . . . , a p } of cardinality p containing exactly all the attractors a i is called the attractor set. Due to the determinism of synchronous Boolean networks, all the attractors are cycles. A cycle is a sequence (x 1 , . . . , x q ) of length q such that ∀j ∈ [ [1, q]], x j+1 = f (x j ) and x q+1 = x 1 : once the system reaches a state x j belonging to a cycle, it successively visits its states x j+1 , . . . , x q , x 1 , . . . , x j for infinity. In the particular case where q = 1, a i is a point attractor. The set B i ⊆ S containing exactly all the x ∈ S from which a i can be reached is called its basin of attraction. With deterministic dynamical systems, the family of sets (B 1 , . . . , B p ) constitutes a partition of S.

Definitions

Some concepts used in this work should be formally defined.

• physiological phenotype: A phenotype which does not impair the life quantity/quality of the organism which exhibits it.

• pathological phenotype: A phenotype which impairs the life quantity/ quality of the organism which exhibits it.

• variant (of a biological network): Given a biological network of interest, a variant is one of its versions, namely the network plus eventually some modifications. It should be noted that this does not exclude the possibility that a variant can be the network of interest as is.

• physiological variant: A variant which produces only physiological phenotypes. It is the biological network of interest as it should be, namely the one of healthy organisms.

• pathological variant: A variant which produces at least one pathological phenotype or which fails to produce at least one physiological phenotype. It is a dysfunctional version of the biological network of interest, namely a version found in ill organisms.

• physiological attractor set: The attractor set A physio of the physiological variant.

• pathological attractor set: The attractor set A patho of the pathological variant.

• physiological Boolean transition function: The Boolean transition function f physio of the physiological variant.

• pathological Boolean transition function: The Boolean transition function f patho of the pathological variant.

• run: An iterative computation of x(k) starting from an x 0 until an a i is reached. It returns w = (x(k 0 ), . . . , x(k end )) where k end depends on when a i is reached, and then on x 0 .

• physiological attractor: An a i such that a i ∈ A physio .

• pathological attractor: An a i such that a i / ∈ A physio .

• modality: The functional perturbation moda i applied on a node v j ∈ V of the network, either activating (moda i = 1) or inactivating (moda i = 0): at each k, moda i overwrites f j (x(k)) making x j (k + 1) = moda i .

• target: A node targ i ∈ V of the network on which a moda i is applied.

• bullet: A couple (c targ , c moda ) where c targ = (targ 1 , . . . , targ r ) is a combination without repetition of targ i and where c moda = (moda 1 , . . . , moda r ) is an arrangement with repetition of moda i , r ∈ [ [1, n]] being the number of targets in the bullet. Here, moda i is intended to be applied on targ i .

• therapeutic bullet: A bullet which makes A patho ⊆ A physio .

• silver bullet: A therapeutic bullet which makes A patho A physio .

• golden bullet: A therapeutic bullet which makes A patho = A physio .

The assumed link between phenotypes and attractors is the reason why attractors are qualified as either physiological or pathological according to the phenotype they produce. This is also the reason why, in this work, target identification aims at manipulating attractor sets of pathological variants.

Steps of the algorithm

The algorithm has two goals: i) finding therapeutic bullets, and ii) classifying them as either golden or silver. A therapeutic bullet makes the pathological variant unable at reaching pathological attractors, that is A patho ⊆ A physio . If such a bullet is applied on a pathological variant, the organism bearing it no longer exhibits the associated pathological phenotypes. However, a therapeutic bullet does not necessarily preserve/restore the physiological attractors. If a therapeutic bullet preserves/restores the physiological attractors, that is if A patho = A physio , then it is a golden one, but if A patho A physio then it is a silver one.

Given a physiological and a pathological variant, that is f physio and f patho , the algorithm follows five steps:

1. with f physio it computes the control attractor set A physio 2. it generates bullets and, for each of them, it performs the three following steps The algorithm is described step by step but can be found as one block of pseudocode in Appendix 6.1 page 33.

Step 1: computing A physio

First of all, A physio has to be computed since it is the control and, as such, determines what is pathological. To do so, runs are performed with f physio and the reached a i are stored in A physio . However, x 0 ∈ S and card S increases exponentially with n. Even for reasonable values of n, card S explodes: more than 1 000 000 possible x 0 for n = 20. One solution ensuring that all the a i are reached is to start a run from each of the possible x 0 , that is from each of the x ∈ S. Practically, this is unfeasible for an arbitrary value of n since the required computational capacity can be too demanding. For example, assuming that a run requires 1 millisecond and that n = 50, performing a run from each of the 2 50 x ∈ S requires nearly 36 000 years.

Given that with deterministic dynamical systems (B 1 , . . . , B p ) is a partition of S, a solution is to select a subset D ⊆ S of a reasonable cardinality containing the x 0 to start from. In this work, D is randomly selected from a uniform distribution. The stumbling block of this solution is that it does not ensure that at least one x 0 per B i is selected and then does not ensure that all the a i are reached. This stumbling block holds only if card D < card S.

Again given that synchronous Boolean networks are deterministic, if a run visits a state already visited during a previous run then its destination, that is the reached attractor, is already found. If so, the run can be stopped and the algorithm can jump to the next one. To implement this, the previous trajectories are stored in a set H, the history, and at each k the algorithm checks if ∃w ∈ H : x(k) ∈ w. If this check is positive then the algorithm jumps to the next run.

To detect the attractors, since with deterministic dynamical systems they are cycles, the algorithm checks at each k if x(k + 1) is an already visited state of the current run, namely if

∃k ′ ∈ [[1, k]] : x(k + 1) = x(k ′ ). If this check is positive then a i = (x(k ′ ), . . . , x(k)).
This step can be written in pseudocode as:

1 prompt card D 2 card D = min(card D, 2 n ) 3 generate D ⊆ S 4 H = {} 5 A physio = {} 6 for x 0 ∈ D do 7 k = 1 8 x(k) = x 0 9
while true do It should be noted that the purpose of this work is not to propose an algorithm for finding Boolean network attractors since advanced algorithms for such tasks are already published [START_REF] Guo | A parallel attractor finding algorithm based on boolean satisfiability for genetic regulatory networks[END_REF][START_REF] Berntenis | Detection of attractors of large boolean networks via exhaustive enumeration of appropriate subspaces of the state space[END_REF][START_REF] Zheng | An efficient algorithm for finding attractors in synchronous boolean networks with biochemical applications[END_REF][START_REF] Dubrova | A sat-based algorithm for finding attractors in synchronous boolean networks[END_REF][START_REF] Ay | Scalable steady state analysis of boolean biological regulatory networks[END_REF]. The purpose is to propose a computational methodology exploiting Boolean network attractors for in silico target identification, a methodology which requires de facto these attractors to be found. This point is discussed in the Conclusion section page 22.

10 if ∃w ∈ H : x(k) ∈ w then 11 break 12 end if 13 x(k + 1) = f physio (x(k)) 14 if ∃k ′ ∈ [[1, k]] : x(k + 1) = x(k ′ ) then 15 A physio = A physio ∪ {(x(k ′ ), . . . ,

Step 2: generating bullets

Bullets are candidate perturbations to apply on the pathological variant to make it unable at reaching pathological attractors and then unable at producing pathological phenotypes. Generating a bullet requires a choice of targ i ∈ V and associated moda i ∈ {0, 1}. In this work, there is no sequencing in target engagement nor in modality application. This means that, given a bullet and during a given run, all the moda i are applied on their corresponding targ i throughout the run. As a consequence, for a given bullet, choosing the same targ i more than once is senseless while it is possible to choose the same moda i for more than one targ i . Therefore, a bullet is a combination c targ without repetition of targ i together with an arrangement c moda with repetition of moda i .

If bullets containing r targets have to be generated then there are n!/(r!•(nr)!) possible c targ and, for each of them, there are 2 r possible c moda . This raises the same computational difficulty than with the state space explosion since there are (n! • 2 r )/(r! • (nr)!) possible bullets. For example, with n = 50 and r = 3 there are more than 150 000 possible bullets. Knowing that the algorithm, as explained below, computes one attractor set per bullet, the computation time becomes practically unfeasible.

To overcome this barrier, the algorithm asks for r as an interval [[r min , r max ]], asks for a maximum number max targ of c targ to generate and asks for a maximum number max moda of c moda to test for each c targ . The algorithm then generates a set C targ of c targ with card C targ ≤ max targ by randomly selecting, from a uniform distribution and without repetition, nodes in the network. In the same way, the algorithm generates a set C moda of c moda with card C moda ≤ max moda by randomly choosing, from a uniform distribution and with repetition, modalities as either activating [START_REF] Knowles | Target selection in drug discovery[END_REF] or inactivating (0). The result is the bullets: per r ∈ [[r min , r max ]], a C targ together with a C moda . As with the state space explosion, the stumbling block of this method is that it does not ensure that all the possible c targ together with all the possible c moda are tested. This stumbling block holds only if

max targ < n!/(r! • (n -r)!) or max moda < 2 r .
This step can be written in pseudocode as:

1 prompt r min , r max , max targ , max moda 2 r max = min(r max , n) 3 golden_set = {} 4 silver_set = {} 5 for r ∈ [[r min , r max ]] do 6 max r targ = min(max targ , n!/(r! • (n -r)!)) 7 max r moda = min(max moda , 2 r ) 8 C targ = {} 9 C moda = {} 10 while card C targ < max r targ do 11 generate c targ / ∈ C targ 12 C targ = C targ ∪ {c targ } 13 end while 14 while card C moda < max r moda do 15 generate c moda / ∈ C moda 16 C moda = C moda ∪ {c moda } 17
end while 18 do steps 3 to 5 19 end for 20 return golden_set, silver_set Line 2 catches the mistake r > n. Lines 3 and 4 create the sets in which the therapeutic bullets found in step 4 are classified as either golden or silver in step 5. Lines 6 and 7 catch the mistake where max targ or max moda is greater than its maximum, which depends on r, hence the creation of max r targ and max r moda to preserve the initially supplied value. Lines 11 and 15 ensure that only new c targ and c moda are generated.

Step 3: computing A patho

Having the control attractor set A physio and a bullet (c targ , c moda ) ∈ C targ × C moda , the algorithm computes the variant attractor set A patho under the effect of (c targ , c moda ) by almost the same way

A physio is computed in step 1. How- ever, f patho is used instead of f physio and (c targ , c moda ) is applied: at each k, f j (x(k)) is overwritten by moda i ∈ c moda , that is x j (k + 1) = moda i , provided that v j = targ i ∈ c targ .
In order to apply all the generated bullets, the algorithm uses two nested f or loops. For each c targ ∈ C targ , it uses successively all the c moda ∈ C moda . For each (c targ , c moda ), the algorithm computes the corresponding A patho and does steps 4 and 5. This step can be written in pseudocode as:

1 for c targ ∈ C targ do 2 for c moda ∈ C moda do 3 H = {} 4 A patho = {} 5 for x 0 ∈ D do 6 k = 1 7 x(k) = x 0 8
while true do end for 30 end for Lines 13-19 are where bullets are applied.

9 if ∃w ∈ H : x(k) ∈ w then 10 break 11 end if 12 x(k + 1) = f patho (x(k)) 13 for targ i ∈ c targ do 14 for v j ∈ V do 15 if v j = targ i then 16 x j (k + 1) = moda i 17 end if 18 end for 19 end for 20 if ∃k ′ ∈ [[1, k]] : x(k + 1) = x(k ′ ) then 21 A patho = A patho ∪ {(x(k ′ ), . . . , x(k))} 22 break 23 end if 24 k = k + 1 25 end while 26 H = H ∪ {(x(1), . . . , x(k))}

Step 4: identifying therapeutic bullets

To identify therapeutic bullets among the generated ones, for each (c targ , c moda ) tested in step 3 and once the corresponding A patho is obtained, the algorithm compares it with A physio to check if A patho ⊆ A physio . This check ensures that the pathological attractors are removed and that if new attractors appear then they are physiological. If this check is positive then the bullet is therapeutic and the algorithm pursues with step 5.

This step can be written in pseudocode as:

1 if A patho ⊆ A physio then 2 do step 5 3 end if 2.3.5
Step 5: assessing therapeutic bullets Therapeutic bullets are qualified as either golden or silver according to their ability at making the pathological variant reaching the physiological attractors. All therapeutic bullets, being golden or silver, remove the pathological attractors without creating new ones, that is A patho ⊆ A physio . However, this does not imply that they preserve/restore the physiological attractors. A golden bullet preserves/restores the physiological attractors: A patho = A physio whereas a silver bullet does not: A patho A physio .

In this setting, golden bullets are perfect therapies whereas silver bullets are not. However, since precious things are rare and just as gold is rarer than silver, finding golden bullets is less likely than finding silver ones. Indeed, given that more constraints are required for a therapeutic bullet to be golden, it is more likely that the found therapeutic bullets are silver, except in one case: card A physio = 1.

Theorem 1. If card A physio = 1 then therapeutic bullets are golden.

Proof.

(therapeutic bullet) ⇒ (A patho ⊆ A physio ) (1) (1) ⇒ (A patho ∈ P(A physio )) (2)
(card A physio = 1) ⇒ (A physio = {a}) (3) (3) ⇒ (P(A physio ) = {∅, {a}}) (4)
((2) ∧ (4)) ⇒ ((A patho = {a}) ∨ (A patho = ∅)) (5) (deterministic dynamical systems) ⇒ (A = ∅)(6) (6) ⇒ (A patho = ∅) (7)
(( 5)

∧ (7)) ⇒ (A patho = {a}) (8) ((3) ∧ (8)) ⇒ (A patho = A physio ) (9)
(9) ⇒ (therapeutic bullet is golden)

Practically, an organism bearing a pathological variant treated with a therapeutic bullet no longer exhibits the associated pathological phenotypes. Moreover, if the therapeutic bullet is golden then the organism exhibits the same phenotypes than its healthy counterpart. However, if the therapeutic bullet is silver then the organism fails to exhibit at least one physiological phenotype. With a silver bullet this is a matter of choice: what is the less detrimental between a silver bullet and no therapeutic bullet at all? This step can be written in pseudocode as:

1 if A patho = A physio then 2 golden_set = golden_set ∪ {(c targ , c moda )} 3 else 4 silver_set = silver_set ∪ {(c targ , c moda )} 5 end if

Example network

To illustrate the algorithm, it is used on a Boolean model of the mammalian cell cycle published by Faure et al [START_REF] Faure | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]. This model is chosen for several reasons: i) a synchronous updating is performed: to date, the algorithm focuses on synchronous Boolean networks, ii) a mammalian biological system is modeled: the closer to human physiology the model is the better it illustrates the intended applications, iii) the cell cycle is a at the heart of cancer: this gives relevancy to the example network, iv) the network comprises ten nodes: easily computable in face of its state space, and v) attractors are already computed: useful to validate the algorithm in finding them. A graphical representation of the example network is shown in Figure 1 page 15. Below are the corresponding Boolean functions where, for the sake of readability, x i stands for x i (k) and x i+ stands for x i (k + 1):

CycD + = CycD Rb + = (¬CycD ∧ ¬CycE ∧ ¬CycA ∧ ¬CycB) ∨ (p27 ∧ ¬CycD ∧ ¬CycB) E2F + = (¬Rb ∧ ¬CycA ∧ ¬CycB) ∨ (p27 ∧ ¬Rb ∧ ¬CycB) CycE + = E2F ∧ ¬Rb CycA + = (E2F ∧ ¬Rb ∧ ¬Cdc20 ∧ ¬(Cdh1 ∧ U bcH10)) ∨(CycA ∧ ¬Rb ∧ ¬Cdc20 ∧ ¬(Cdh1 ∧ U bcH10)) p27 + = (¬CycD ∧ ¬CycE ∧ ¬CycA ∧ ¬CycB) ∨(p27 ∧ ¬(CycE ∧ CycA) ∧ ¬CycB ∧ ¬CycD) Cdc20 + = CycB Cdh1 + = (¬CycA ∧ ¬CycB) ∨ Cdc20 ∨ (p27 ∧ ¬CycB) U bcH10 + = ¬Cdh1 ∨ (Cdh1 ∧ U bcH10 ∧ (Cdc20 ∨ CycA ∨ CycB)) CycB + = ¬Cdc20 ∧ ¬Cdh1
Having the example network, two variants are needed: the physiological one and the pathological one. The physiological variant is the network as is while the pathological variant is the network plus a constitutive activation/inactivation of at least one of its nodes. For simplicity, and given the relatively small number of entities, only one is chosen: the retinoblastoma protein Rb for which a constitutive inactivation is applied. To implement this, the corresponding f i becomes:

Rb(k + 1) = 0 in f patho .
Rb is chosen because its inactivation occurs in many cancers [START_REF] Charles | The rb and p53 pathways in cancer[END_REF]. Therefore, a network bearing a constitutive inactivation of it should be a relevant example of a pathological variant.

Case study

To illustrate the intended usage of the proposed methodology, the algorithm is used on a Boolean model of the Fanconi Anemia/Breast Cancer (FA/BRCA) pathway published by Rodriguez et al [START_REF] Rodriguez | A boolean network model of the fa/brca pathway[END_REF]. This model is chosen for several reasons: i) two pathological conditions are studied: required for a case study of an in silico target identification, ii) the physiological and pathological variants are clearly described: required by the algorithm, iii) it is nearly three times bigger than the example network: representative of a more comprehensive biological model while remaining computationally tractable, iv) synchronous updating is used: to date, the algorithm focuses on synchronous Boolean networks, and v) attractors are already interpreted in terms of phenotypes.

The FA/BRCA pathway is dedicated to DNA repair, more precisely to interstrand cross-link (ICL) removal. As expected with any DNA repair impairment, individuals suffering from FA/BRCA pathway malfunction are subjected to increased risk of cancer, such as in Fanconi anemia, a rare genetic disorder causing bone marrow failure, congenital abnormalities and increased risk [START_REF] Faure | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]. CDKs (cyclin-dependent kinases) are the catalytic partners of cyclins and, in this model, are not explicitly shown since the activity of CDK-cyclin complexes essentially depends on cyclins. Furthermore, the inhibition of E2F by Rb is modeled by opposing Rb to the effects of E2F on its targets. The same applies to the inhibition of CycE and CycA by p27. For a complete description of the model, see [START_REF] Faure | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]. CycD: CDK4/6-cyclin D complex, input of the model, initiates the cell cycle, activated by positive signals such as growth factors; CycE: CDK2-cyclin E complex; CycA: CDK2-cyclin A complex; CycB: CDK1-cyclin B complex; Rb: retinoblastoma protein, a tumor suppressor; E2F: a family of transcription factors divided into activator and repressor members, in this model E2F represents the activator members; p27: p27/Kip1, a CKI (CDK inhibitor); Cdc20: an APC (Anaphase Promoting Complex, an E3 ubiquitin ligase) activator; Cdh1: an APC activator; UbcH10: an E2 ubiquitin conjugating enzyme.

of cancer [START_REF] Robert | Susceptibility pathways in fanconi's anemia and breast cancer[END_REF][START_REF] Auerbach | Fanconi anemia and its diagnosis[END_REF][START_REF] Winter | The genetic and molecular basis of fanconi anemia[END_REF]. Rodriguez et al propose a Boolean model comprising the FA/BRCA pathway and three types of DNA damages commonly observed in Fanconi anemia, namely ICLs, double-strand breaks (DSBs) and DNA adducts (ADDs). It should be noted that the ICL repair process creates DSBs and ADDs before removing them, thus leaving an undamaged DNA ready for the cell cycle. For a complete description of the model, see [START_REF] Rodriguez | A boolean network model of the fa/brca pathway[END_REF]. Below are the corresponding Boolean functions where, for the sake of readability, x i stands for x i (k) and x i+ stands for x i (k + 1):

ICL + = ICL ∧ ¬DSB F AN CM + = ICL ∧ ¬CHKREC F Acore + = F AN CM ∧ (AT R ∨ AT M ) ∧ ¬CHKREC F AN CD2I + = F Acore ∧ ((AT M ∨ AT R) ∨ (H2AX ∧ DSB)) ∧ ¬U SP 1 M U S81 + = ICL F AN CJBRCA1 + = (ICL ∨ ssDN ARP A) ∧ (AT M ∨ AT R) XP F + = (M U S81 ∧ p53 ∧ ¬(F Acore ∧ F AN CD2I ∧ F AN 1)) ∨(M U S81 ∧ ¬F AN CM ) F AN 1 + = M U S81 ∧ F AN CD2I ADD + = (ADD ∨ (M U S81 ∧ (F AN 1 ∨ XP F ))) ∧ ¬P CN AT LS DSB + = (DSB ∨ F AN 1 ∨ XP F ) ∧ ¬(N HEJ ∨ HRR) P CN AT LS + = (ADD ∨ (ADD ∧ F Acore)) ∧ ¬(U SP 1 ∨ F AN 1) M RN + = DSB ∧ AT M ∧ ¬((KU ∧ F AN CD2I) ∨ RAD51 ∨ CHKREC) BRCA1 + = DSB ∧ (AT M ∨ CHK2 ∨ AT R) ∧ ¬CHKREC ssDN ARP A + = DSB ∧ ((F AN CD2I ∧ F AN CJBRCA1) ∨ M RN ) ∧ ¬(RAD51 ∨ KU ) F AN CD1N + = (F AN CD2I ∧ ssDN ARP A ∧ ¬CHKREC) ∨(ssDN ARP A ∧ BRCA1) RAD51 + = ssDN ARP A ∧ F AN CD1N ∧ ¬CHKREC HRR + = DSB ∧ RAD51 ∧ F AN CD1N ∧ BRCA1 ∧ ¬CHKREC U SP 1 + = ((F AN CD1N ∧ F AN CD2I) ∨ P CN AT LS) ∧ ¬F AN CM KU + = DSB ∧ ¬(M RN ∨ F AN CD2I ∨ CHKREC) DN AP K + = (DSB ∧ KU ) ∧ ¬CHKREC N HEJ + = ((DSB ∧ DN AP K ∧ KU ) ∧ ¬(AT M ∧ AT R)) ∨(¬((F AN CJBRCA1 ∧ ssDN ARP A) ∨ CHKREC) ∧DSB ∧ DN AP K ∧ XP F ) AT R + = (ssDN ARP A ∨ F AN CM ∨ AT M ) ∧ ¬CHKREC AT M + = (AT R ∨ DSB) ∧ ¬CHKREC p53 + = (((AT M ∧ CHK2) ∨ (AT R ∧ CHK1)) ∨ DN AP K) ∧ ¬CHKREC CHK1 + = (AT M ∨ AT R ∨ DN AP K) ∧ ¬CHKREC CHK2 + = (AT M ∨ AT R ∨ DN AP K) ∧ ¬CHKREC H2AX + = DSB ∧ (AT M ∨ AT R ∨ DN AP K) ∧ ¬CHKREC CHKREC + = ((P CN AT LS ∨ N HEJ ∨ HRR) ∧ ¬DSB) ∨((¬ADD) ∧ (¬ICL) ∧ (¬DSB) ∧ ¬CHKREC)
The physiological variant is the FA/BRCA pathway model as is. To it, Rodriguez et al propose two pathological variants, here called patho1 and patho2, modeling two mutations involving genes of the FA/BRCA pathway. These mutations are observed in patients suffering from Fanconi anemia [START_REF] Neveling | Genotype-phenotype correlations in fanconi anemia[END_REF]. The first one involves the FANCA gene, corresponding to the F Acore variable, and the second one involves the FANCD1/BRCA2 or FANCN/PALB2 gene, corresponding to the F AN CD1N variable. These mutations are of loss-of-function kind: to simulate them, the corresponding f i become F Acore(k + 1) = 0 for FANCA gene null mutation in f patho1 and F AN CD1N (k + 1) = 0 for FANCD1/BRCA2 or FANCN/PALB2 gene null mutation in f patho2 .

Implementation

The algorithm is implemented in Fortran compiled with GFortran 1 . The code is available on GitHub 2 at https://github.com/arnaudporet/kali-targ.

Results

In this section, results produced with the algorithm on the example network are exposed to illustrate how it works. Next, results produced with the algorithm on the case study are exposed to illustrate its intended applications for target identification.

Results of step 1

Owing to the relatively small size of the example network, card D is set to card S = 1024. Since card D = card S, all the attractors are found. Attractors are presented as matrices where, for an attractor of length q, lines correspond to the x i (k), k ∈ [ [1, q]], and columns to x(k). The algorithm returns the following attractors:

a 1 = CycD 1 1 1 1 1 1 1 Rb 0 0 0 0 0 0 0 E2F 0 1 1 1 0 0 0 CycE 0 0 1 1 1 0 0 CycA 0 0 0 1 1 1 1 p27 0 0 0 0 0 0 0 Cdc20 1 0 0 0 0 0 1 Cdh1 1 1 1 1 0 0 0 U bcH10 1 1 0 0 0 1 1 CycB 0 0 0 0 0 1 1 a 2 = CycD 0 Rb 1 E2F 0 CycE 0 CycA 0 p27 1 Cdc20 0 Cdh1 1 U bcH10 0 CycB 0 1 http://

Results of steps 2 to 5

Results of steps 2 to 5 are grouped since only the therapeutic bullets found in step 4 and classified in step 5 are returned. The algorithm is launched with r min = 1 and r max = 2. Due to the relatively small size of the example network, max targ and max moda are set to their maximum, namely max targ = 45 and max moda = 4. Consequently, all the possible bullets made of 1 to 2 targets are tested. The algorithm returns the following therapeutic bullets: +CycD silver +CycD -p27 silver -CycD +Rb silver +CycD -Rb silver where + means therapeutic activation andmeans therapeutic inactivation. It should be noted that no golden bullets are found, an unsurprising result since they are rarer than silver ones.

Given these results, therapeutic activation of Rb, which is inactivated in the pathological variant, is not enough to remove the pathological attractors. Indeed, as seen in the third bullet, therapeutic activation of Rb must be accompanied by therapeutic inactivation of CycD. To better illustrate what is performed to obtain these therapeutic bullets, below is A patho without any bullet:

a 3 = CycD 0 0 0 0 0 0 0 0 Rb 0 0 0 0 0 0 0 0 E2F 1 1 1 1 0 0 0 0 CycE 0 1 1 1 1 0 0 0 CycA 0 0 1 1 1 1 1 0 p27 1 1 1 0 0 0 0 0 Cdc20 0 0 0 0 0 0 1 1 Cdh1 1 1 1 1 0 0 0 1 U bcH10 1 0 0 0 0 1 1 1 CycB 0 0 0 0 0 1 1 0 a 4 = CycD 1 1 1 1 1 1 1 Rb 0 0 0 0 0 0 0 E2F 1 1 1 0 0 0 0 CycE 0 1 1 1 0 0 0 CycA 0 0 1 1 1 1 0 p27 0 0 0 0 0 0 0 Cdc20 0 0 0 0 0 1 1 Cdh1 1 1 1 0 0 0 1 U bcH10 1 0 0 0 1 1 1 CycB 0 0 0 0 1 1 0
each of these two attractors attracting 50% of the x ∈ S under f patho . It should be noted that a 4 = a 1 ∈ A physio : a 4 is a physiological attractor which also belongs to A patho . Indeed, it is possible that the pathological variant exhibits physiological attractors: A patho is not the set containing exactly all the pathological attractors, it is the attractor set of the pathological variant, so A physio ∩ A patho = ∅ is possible. However, a 3 / ∈ A physio : it is a pathological attractor and is what a therapeutic bullet, being golden or silver, is intended to remove.

Again to better illustrate what is performed to obtain these therapeutic bullets, below is A patho under the third bullet:

CycD 0 Rb 1 E2F 0 CycE 0 CycA 0 p27 1 Cdc20 0 Cdh1 1 U bcH10 0 CycB 0 which is a 2 .
As expected for a therapeutic bullet, the pathological attractor a 3 is removed. However, the physiological attractor a 1 is not restored: the third therapeutic bullet is silver. Consequently, with this therapeutic bullet no cell cycle occurs and the only reachable phenotype is quiescence. While disabling the cell cycle of cancer cells is beneficial, disabling the cell cycle of healthy cells is not. As mentioned above, with silver bullets this is a matter of choice.

Results of the case study

With the case study, card S = 268 435 456: computing attractors from all the x ∈ S becomes too demanding. Indeed, it should be recalled that the algorithm computes one attractor set per bullet, namely A patho under the tested bullet. Consequently, card D is set to a more reasonable value: card D = 10 000. Despite that card D < card S, it seems sufficient for the algorithm to find all the attractors, just as Rodriguez et al whose the computation covers the whole state space. Below are the computed attractors: where

• A physio = {a 1 } • A patho1 = {a 1 } • A patho2 = {a 1 , a 2 },
a 1 = ICL 0 0 F AN CM 0 0 F Acore 0 0 F AN CD2I 0 0 M U S81 0 0 F AN CJBRCA1 0 0 XP F 0 0 F AN 1 0 0 ADD 0 0 DSB 0 0 P CN AT LS 0 0 M RN 0 0 BRCA1 0 0 ssDN ARP A 0 0 F AN CD1N 0 0 RAD51 0 0 HRR 0 0 U SP 1 0 0 KU 0 0 DN AP K 0 0 N HEJ 0 0 AT R 0 0 AT M 0 0 p53 0 0 CHK1 0 0 CHK2 0 0 H2AX 0 0 CHKREC 0 1 a 2 = ICL 0 F AN CM 0 F Acore 0 F AN CD2I 0 M U S81 0 F AN CJBRCA1 1 XP F 0 F AN 1 0 ADD 0 DSB 1 P CN AT LS 0 M RN 1 BRCA1 1 ssDN ARP A 1 F AN CD1N 0 RAD51 0 HRR 0 U SP 1 0 KU 0 DN AP K 0 N HEJ 0 AT R 1 AT M 1 p53 1 CHK1 1 CHK2 1 H2AX 1 CHKREC 0
and their biological interpretation:

• a 1 : cell cycle progression

• a 2 : cell cycle arrest

In physiological conditions, in case of a damaged DNA, cells repair it before performing the cell cycle, or die if repair fails. Such checkpoints enable cells to ensure genomic integrity by preventing damaged DNA to be replicated and then propagated [START_REF] Bartek | Dna damage checkpoints: from initiation to recovery or adaptation[END_REF][START_REF] Ishikawa | Dna damagedependent cell cycle checkpoints and genomic stability[END_REF]. Otherwise, genetic instability may appears, potentially leading to cancer [START_REF] Nakanishi | Genetic instability in cancer cells by impaired cell cycle checkpoints[END_REF]. The results show that the physiological variant is able to ensure genomic integrity since its unique attractor is a 1 where ICL = DSB = ADD = 0: DNA damages are repaired, if any, and the cell cycle can safely occur. Interestingly, the same physiological phenotype is computed for patho1 where A patho1 = A physio . This suggests that cells bearing FANCA gene null mutation are nonetheless able to repair DNA. With patho2, a pathological attractor appears: a 2 , where DSB = 1. This suggests that cells bearing FANCD1/BRCA2 or FANCN/PALB2 gene null mutation are unable to repair DSBs, explaining why a 2 corresponds to cell cycle arrest: DNA remains damaged. It should be noted that a 1 ∈ A patho2 , suggesting that from some x 0 , that is under some conditions, such cells could be able to repair DNA. However, a 1 attracts only 29.5% of the x ∈ D under f patho2 , indicating that the pathological phenotype associated with a 2 is the most likely.

Altogether, according to the computed attractors and their phenotypic interpretation, and limited to the scope studied by the model of Rodriguez et al, FANCA gene null mutation may not induce pathological phenotypes. However, with FANCD1/BRCA2 or FANCN/PALB2 gene null mutation, two phenotypes are predicted: a physiological one and a pathological one, the latter being the most likely. Therefore, the algorithm has to operate on patho2 to find bullets able to remove the pathological attractor a 2 . By comprehensively testing all the bullets made of 1 to 3 targets, the algorithm returns the following results: number of all possible bullets number of therapeutic bullets r = 1 In this case study, DNA damages such as ICLs and DSBs are the pathological events. Unsurprisingly, the algorithm suggests them to be targeted: this is a logical consequence. However, DNA damages are not biomolecules in themselves and directly targeting them by means of drugs appears senseless. What is relevant are the biomolecules of the FA/BRCA pathway suggested as therapeutic targets. Interestingly, ATM dominates all the other candidates, predicting it to be a pivotal therapeutic target for the patho2 condition, namely the FA/BRCA pathway bearing FANCD1/BRCA2 or FANCN/PALB2 gene null mutation, as observed in Fanconi anemia.

Conclusion

Under the assumption that attractors of dynamical systems and phenotypes of biological networks are linked when the former models the latter, the results show that the algorithm succeeds in performing the proposed in silico target identification. It returns therapeutic bullets for a pathological variant of the mammalian cell cycle relevant in cancer and for a pathological variant modeling Fanconi anemia. Consequently, the algorithm can be used on other synchronous Boolean models of biological networks involved in diseases for in silico target identification. It is intended to be of use in the early steps of target identification by providing an efficient way to identify candidate targets prior to costly wet lab experiments. However, both the physiological and pathological variants have to be known. This can constitute a limit of the proposed methodology since not all the pathophysiologies are known. Target identification, whether performed in silico or not, is a step belonging to a wider process: drug discovery. Having demonstrated a potential target in silico, or even in vitro, is far from having a medication. Further work and many years are necessary before obtaining a drug which is effective in vivo. For example, and among other characteristics, such a drug has to be absorbed by the organism, has to reach its target and has to be non-toxic at therapeutic dosages. Furthermore, as with any in silico evidence, it should be validated through wet lab experiments: there is a bridge to cross between theory and practice. Indeed, mathematical models approximate reality without reproducing it and theory must meet practice. For example, targeting ATM should restore a physiological running of the FA/BRCA pathway bearing FANCD1/BRCA2 or FANCN/PALB2 gene null mutation. However, if ATM operates in other pathways, targeting it may disturb them, thus potentially creating de novo non-physiological conditions. Nevertheless, it is expected that the algorithm is of interest for target identification, notably by exploiting the inexpensiveness and predictive power of computational approaches to optimize the efficiency of costly wet lab experiments.. While finding Boolean network attractors of biological networks is not the purpose of this work, it is a necessary step which is in itself a challenging field of computational biology. Therefore, incorporating advances made in this field could be an interesting improvement. Another possible improvement could be to extend the algorithm for asynchronous Boolean networks since such models are likely to more accurately describe the dynamics of biological systems [START_REF] Zhu | Asynchronous stochastic boolean networks as gene network models[END_REF][START_REF] Liang | Stochastic boolean networks: an efficient approach to modeling gene regulatory networks[END_REF]. Indeed, in biological systems, events may be subjected to stochasticity, may not occur simultaneously or may not belong to the same time scale, three points that a synchronous updating scheme does not take into account.

Yet another possible improvement could be to use a finer logic, such as multivalued logic. One of the main limitations of Boolean models is that variables can take only two values. In reality, things are not necessarily binary and variables should be able to take more values. Multivalued logic enables it in a discrete manner where variables can take a finite number of values between 0 (false) and 1 (true). For example, one can state that Rb is partly impaired rather than totally. Such a statement is not implementable with Boolean models but is with multivalued ones such as, for example, a three-valued logic where true = 1, moderate = 0.5 and f alse = 0.

Finally, considering the basin of attraction of the pathological attractors could be an interesting extension of the criterion for selecting therapeutic bullets. In that case, the therapeutic potential of bullets could be assessed by estimating their ability at reducing the basin of the pathological attractors, as performed by Fumia et al with their Boolean model of cancer pathways [START_REF] Herman | Boolean network model for cancer pathways: predicting carcinogenesis and targeted therapy outcomes[END_REF]. Such a criterion enables to consider the particular case where pathological attractors are removed, that is where pathological basins are reduced to the empty set, but also the other cases where pathological basins are not necessarily reduced to the empty set. Such a less restrictive selection of therapeutic bullets would enable to consider more targeting strategies for counteracting diseases.

Additional improvements

First of all, some additional definitions should be stated:

• physiological state space: The state space S physio of the physiological variant.

• pathological state space: The state space S patho of the pathological variant.

• testing state space: The state space S test of the pathological variant under the effect of a bullet.

• physiological basin: The basin of attraction B physio,i of a physiological attractor a physio,i .

• pathological basin: The basin of attraction B patho,i of a pathological attractor a patho,i .

• n-bullet: A bullet made of n targets.

Among the possible improvements mentioned in the Conclusion section page 22, two are done: extending the algorithm for multivalued logic and considering pathological basins for selecting therapeutic bullets. Therefore, variables of multivalued logic can model more than only two possible levels, enabling models to be more realistic than those based on Boolean logic.

Multivalued logic

Methods

Boolean logic can be seen as a particular case of multivalued logic: it is a bivalued logic where variables take their value in {0, 1}. While Boolean operators work well in this case, multivalued logic requires suitable logical operators to be introduced. One solution is to use a mathematical formulation of the Boolean operators which also works with any multivalued logic, just as the Zadeh operators. These logical operators are a mathematical generalization of the Boolean ones proposed for fuzzy logic by its pioneer Lotfi Zadeh. Their mathematical formulation is as follow:

AN D(x, y) = min(x, y) OR(x, y) = max(x, y) N OT (x) = 1 -x
With a h-valued logic, card S = h n . If h = 2 then this is the Boolean case, where card S already raises computational difficulties. With an arbitrary h > 2, card S raises even more computational difficulties. The same applies to the testable bullets since there are h r possible c moda and then (n!•h r )/(r!•(n-r)!) possible bullets. To illustrate how the algorithm works with a multivalued logic without overloading it, a 3-valued logic is used with {0, 0.5, 1} as domain of value: x i (k) ∈ {0, 0.5, 1}. 0 and 1 have the same meaning as in Boolean logic, namely false and true respectively. 0.5 is an intermediate truth degree which can be seen as an intermediate level of activity or abundance, depending on what is modeled. Consequently, S = {0, 0.5, 1} n implying x 0 , x(k) ∈ {0, 0.5, 1} n , D ⊆ {0, 0.5, 1} n and moda i ∈ {0, 0.5, 1}. Moreover, the Boolean operators of the f i are replaced by the Zadeh operators. This results in the following minor changes in the pseudocode of the algorithm described in Appendix 6.1 page 33:

line Boolean logic h-valued logic 2 card D = min(card D, 2 n ) card D = min(card D, h n ) 29 max r moda = min(max moda , 2 r ) max r moda = min(max moda , h r )
How the algorithm works with this 3-valued logic is illustrated with the example network, whose the logical functions become:

CycD + = CycD Rb + = max(min(1 -CycD, 1 -CycE, 1 -CycA, 1 -CycB), min(p27, 1 -CycD, 1 -CycB)) E2F + = max(min(1 -Rb, 1 -CycA, 1 -CycB), min(p27, 1 -Rb, 1 -CycB)) CycE + = min(E2F, 1 -Rb) CycA + = max(min(E2F, 1 -Rb, 1 -Cdc20, 1 -min(Cdh1, U bcH10)), min(CycA, 1 -Rb, 1 -Cdc20, 1 -min(Cdh1, U bcH10))) p27 + = max(min(1 -CycD, 1 -CycE, 1 -CycA, 1 -CycB), min(p27, 1 -min(CycE, CycA), 1 -CycB, 1 -CycD)) Cdc20 + = CycB Cdh1 + = max(min(1 -CycA, 1 -CycB), Cdc20, min(p27, 1 -CycB)) U bcH10 + = max(1 -Cdh1, min(Cdh1, U bcH10, max(Cdc20, CycA, CycB))) CycB + = min(1 -Cdc20, 1 -Cdh1)
which is f physio . For f patho , owing to this 3-valued logic, a constitutive but partial inactivation of Rb is simulated. The corresponding f i becomes:

Rb + = 0.5 in f patho .

Results

With the example network modeled by this 3-valued logic, card S = 59 049, which remains computationally tractable. Therefore, card D = card S: all the attractors are found. With the physiological variant, the algorithm returns:

A physio = {a physio1 , a physio2 ,
a physio2 = CycD 0 Rb 1 E2F 0 CycE 0 CycA 0 p27 1 Cdc20 0 Cdh1 1 U bcH10 0 CycB 0 a physio3 = CycD
a physio6 = CycD 1 1 1 1 1 1 Rb 0 0 0 0 0 0 E2F 0 1 1 1 0 0 CycE 0 0 1 1 1 0 CycA 0 0 0 1 1 1 p27 0 0 0 0 0 0 Cdc20 1 0 0 0 0 0 Cdh1 1 1 1 1 0 0 U bcH10 1 1 0 0 0 1 CycB 0 0 0 0 0 1
and their corresponding basin of attraction:

a i B i (in % of card S physio ) a physio1 9.9% a physio2 20.1% a physio3
33.3% a physio4

24.5% a physio5

3.4% a physio6

8.8%

It should be noted that a physio2 and a physio6 are the two physiological attractors found in the Boolean case. Indeed, since {0, 1} ⊂ {0, 0.5, 1} and since the Zadeh operators also work with Boolean logic, Boolean logic is included in this three-valued logic. This means that results obtainable with the former are also obtainable with the latter. With the pathological variant, where Rb is constitutively but partially inactivated, the algorithm returns:

A patho = {a physio1 , a physio3 , a patho1 }
where

a patho1 = CycD 1 Rb 0.5 E2F 0.5 CycE 0.5 CycA 0.5 p27 0 Cdc20 0.5 Cdh1 0.5 U bcH10 0.5 CycB 0.5
ans their corresponding basin of attraction:

a i B i (in % of card S patho ) a physio1
33.3% a physio3

33.3% a patho1

33.3%

Only a physio1 and a physio3 remain, while a patho1 appears and is what therapeutic bullets have to remove from S test .

As in the Boolean case, the algorithm is launched with r min = 1 and r max = 2. max targ and max moda are set to their maximum, namely max targ = 45 and max moda = 9: all the 1, 2-bullets are tested. The algorithm returns the following therapeutic bullets:

CycD[0] silver CycD[0.5] silver CycD[0] Rb[0.5] silver CycD[0.5] Rb[0.5] silver CycD[1] Rb[0] silver CycD[0] E2F [0.5] silver CycD[0.5] E2F [0.5] silver CycD[0] CycE[0.5] silver CycD[0.5] CycE[0.5] silver CycD[0] CycA[0.5] silver CycD[0.5] CycA[0.5] silver CycD[0] p27[0.5] silver CycD[0.5] p27[0.5] silver CycD[0] Cdc20[0.5] silver CycD[0.5] Cdc20[0.5] silver CycD[0] Cdh1[0.5] silver CycD[0.5] Cdh1[0.5] silver CycD[0] U bcH10[0.5] silver CycD[0.5] U bcH10[0.5] silver CycD[0] CycB[0.5] silver CycD[0.5] CycB[0.5] silver
where X[y] means that the node X ∈ V has to be set to the value y ∈ {0, 0.5, 1}.

For example, the third therapeutic bullet is made of the targets CycD and Rb whose the value has to be set to 0 and 0.5 respectively. As in the Boolean case, it should be noted that no golden bullets are found, an unsurprising result since they are rarer than silver ones.

Conclusion

The algorithm is now extended for multivalued logic, which includes the Boolean one. This means that the previous strictly Boolean version of the algorithm is included in this new one. Moreover, allowing variables to take an arbitrary finite number of values should enable to more accurately model biological processes and produce more fine-tuned therapeutic bullets. However, this accuracy and fine-tuning are at the cost of an increased computational requirement. Indeed, in this work, the computational requirement essentially depends on the cardinality of the state space, which itself depends on the size of the model and the used multivalued logic. Therefore, the size of the model and the used multivalued logic should be balanced: the smaller the model is, the more variables should be finely valued. For example, for a fine therapeutic investigation, the model should only contain the essential and specific pieces of the pathophysiology of interest, modeled by a finely valued logic. On the other hand, for a gross therapeutic investigation, an exhaustive model could be used but modeled by a coarse-grained logic, such as the Boolean one. Finally, it should be noted that the ultimate multivalued logic is the infinitely valued one, which is fuzzy logic. With fuzzy logic, the whole [0; 1] ⊂ R is used to valuate variables, which should bring the best accuracy for the qualitative modeling formalism [START_REF] Poret | A logicbased modeling derived from boolean networks: adding fuzzy logic and edge tuning[END_REF].

Therapeutic bullet assessment

Introduction

Till now, the algorithm requires therapeutic bullets to remove all the pathological attractors from the pathological state space, so that the pathological variant no longer exhibits pathological phenotypes. This criterion for selecting therapeutic bullets can appear somewhat drastic since it is all or nothing. A less strict criterion should enable to consider more targeting strategies, and then more possibilities for counteracting diseases. Certainly, a less restrictive criterion could bring less "powerful" therapeutic bullets, but being too demanding potentially leads to no results and loss of nonetheless interesting findings. The therapeutic potential of bullets could be assessed by estimating their ability at reducing the cardinality of the pathological basins. This is a more permissive criterion since therapeutic bullets no longer have to necessarily remove the pathological attractors. Reducing the cardinality of a pathological basin renders the corresponding pathological attractor less reachable, and then the associated pathological phenotype less likely. This new criterion includes the previous one: removing an attractor means reducing its basin of attraction to the empty set. Therefore, therapeutic bullets obtainable with the previous criterion are also obtainable with this new one.

Methods

To implement this new criterion for selecting therapeutic bullets, the algorithm considers a bullet as therapeutic if it increases card B physio,i in S test without creating de novo attractors. Since the attractors are either physiological or pathological, increasing card B physio,i is equivalent to decreasing card B patho,i . The goal of this new criterion is to increase the physiological part of S test , which is equivalent to decreasing its pathological part. Consequently, a pathological variant treated by such a therapeutic bullet tends to, but not necessarily reaches, an overall physiological behavior. However, as with the previous criterion, it does not ensure that the a physio,i are preserved/restored. A fortiori, it does not ensure that the B physio,i in S test are as in S physio . This means that it does not ensure that the reachability of the a physio,i is preserved/ restored. Nevertheless, as with the previous criterion, this is a matter of choice between a therapeutic bullet or not. To assist this choice and better visualize the effects of therapeutic bullets, the card B physio,i and card B patho,i in S test are computed.

Implementing this new criterion for selecting therapeutic bullets is a major change. Therefore, the pseudocode of the algorithm presented in Appendix 6.1 page 33 is rewritten and structured into three modules:

• the compute_A function, which computes A physio or A patho , depending on which of the f physio or f patho is passed

• the compute_cover function, which for two attractor sets A 1 and A 2 computes the covering of S 2 by B 1,i , expressed in percents of card S 2

• the compute_T function, which computes a set T of therapeutic bullets Below is the corresponding pseudocode: For A physio and A patho , which are computed without bullet, the empty bullet ((), ()) has to be passed. The a i are represented as structures composed of two fields: a i .seq, which is the sequence of a i (line 15), and a i .f req, which is the corresponding card B i , expressed in percents of card D. To compute a i .f req, the algorithm counts the number of times a i is reached (line 19 if this is the first time a i is reached, line 17 otherwise) and then, once all the x 0 ∈ D are computed, translates a i .f req in percents of card D (line 28).

function y = compute_cover(A 1 , A 2 ) 1 cover = 0 2 for a 1 ∈ A 1 do 3 if ∃a 2 ∈ A 2 : a 1 .seq = a 2 .seq then 4 cover = cover + a 2 .f req 5
end if 6 end for 7 return cover end function If a 1 also belongs to A 2 (line 3) then the cardinality of its basin in S 2 is used to compute the covering of S 2 by B 1,i (line 4). function T = compute_T (f physio , f patho , r min , r max , max targ , max moda , max D , h, V ) end for 33 end for 34 return T end function max D is the desired card D and h is the cardinality of the domain of value, which depends on the used multivalued logic. A physio and A patho are computed without bullet, so the empty bullet ((), ()) is passed to compute_A (lines 7 and 8). cover patho is the covering of S patho by B physio,i (line 10) and cover test is the covering of S test by B physio,i (line 26). A test is the pathological attractor set under the effect of the tested bullet (line 24). A therapeutic bullet has to avoid the appearance of de novo attractors (line 25) and has to increase the covering of S test by B physio,i (line 27).

1 n = card V 2 D = {} 3 while card D < max D do 4 generate x 0 / ∈ D 5 D = D ∪ {x 0 } 6 end while 7 A physio = compute_A(f physio , (), (), D, V ) 8 A patho = compute_A(f patho , (), (), D, V ) 9 T = {} 10 cover patho = compute_cover(A physio , A patho ) 11 for r ∈ [[r min ,

Results

This new criterion for selecting therapeutic bullets is illustrated on the case study modeled by Boolean logic: h = 2. Since patho1 has the same attractor set than the physiological variant, only patho2 is computed. As previously, wholly computing S is too demanding. Therefore, D is intended to have a reasonable cardinality: max D = 100 000. All the 1, 2-bullets are tested: r min = 1, r max = 2, max targ = 378 and max moda = 4. However, their therapeutic potential is no longer expressed as golden or silver but by their gain. It is displayed as follow: x% → y% where card B physio,i = x% in S patho and y% in S test . Consequently, in order to increase the physiological part of S test , a therapeutic bullet has to make y > x. The card B physio,i and card B patho,i in S test are also computed and expressed in percents of card S test . The algorithm returns 59 therapeutic bullets whose the list can be found in Appendix 6.3 page 39.

A therapeutic bullet as defined by the previous criterion, that is which removes all the a patho,i from S test , makes de facto card B physio,i = 100% in S test . As already mentioned, the previous criterion is included in this new one: therapeutic bullets obtainable with the former are also obtainable with the latter. This can be checked by noting that the 1, 2-therapeutic bullets found with the previous criterion are also found with this new one.

With this case study, A physio = {a physio1 }, so B physio,i = B physio1 . Therefore, in this particular case where card A physio = 1, therapeutic bullets have to increase card B physio1 in S test . It should be recalled that card B physio1 = 29.4% in S patho , so therapeutic bullets have to make card B physio1 > 29.4% in S test .

For example, below are the computed -AT M is a therapeutic bullet also found with the previous criterion since it removes all the a patho,i , namely a patho1 , from S test . However, the other four therapeutic bullets are only obtainable with this new criterion since they do not remove a patho1 from S test . Nevertheless, as therapeutic bullets, they increase card B physio1 in S test . This highlight the ability of this new criterion to unravel more therapeutic bullets of varying therapeutic potential, thus opening the way for more targeting strategies of varying theoretical efficacy. Of course, therapeutic bullets of poor potential are also unraveled, such as -F AN CD2I which only increases card B physio1 from 29.4% in S patho to 30.4% in S test . However, in silico tools should not restrict their predictions to only those exhibiting a high theoretical potency since predicted does not necessarily mean true. Indeed, a prediction of apparently poor interest can reveal itself of great interest in practice, and vice versa.

Conclusion

The algorithm now uses a new criterion for selecting therapeutic bullets which brings a wider range of targeting strategies intended to push pathological behaviors toward physiological ones with varying predicted efficacy. Moreover, no information is lost from the previous criterion since results obtainable with the previous one are also obtainable with this new one. This new criterion is based on a more permissive assumption stating that reducing the reachability of pathological attractors is therapeutic. For an in silico tool such as this algorithm, a more permissive assumption is important since theoretical findings have to outlive the bottleneck separating prediction to reality. With a too strict criterion, the risk of highlighting too few candidate targets or to miss some interesting ones is too hight. Indeed, results predicted in silico have to be validated in vitro and/or in vivo. Therefore, requiring only perfect predictions such as therapeutic bullets removing all the pathological attractors could left insufficient results after validation. All the more so that a prediction of apparently poor interest could reveal itself as an insight of great interest and vice versa.

This new criterion for selecting therapeutic bullets also brings a finer assessment of their potential since all the percentages between card B physio,i in S patho and 100% are considered. With the previous criterion, the only therapeutic potential is card B physio,i = 100% in S test , thus reducing the assessment to therapeutic or not. However, things are not necessarily black or white but rather a continuum of gray nuances, so the assessment of therapeutic potentials should be nuanced too, just as it is now.

Appendices

Appendix 1

The algorithm in one block of pseudocode.

1 prompt card D 2 card D = min(card D, 2 n ) 3 generate D ⊆ S 4 H = {} 5 A physio = {} 6 for x 0 ∈ D do 7 k = 1 8 x(k) = x 0 9 while true do 10 if ∃w ∈ H : x(k) ∈ w then 11 break 12 end if 13 x(k + 1) = f physio (x(k)) 14 if ∃k ′ ∈ [[1, k]] : x(k + 1) = x(k ′ ) then 15 A physio = A physio ∪ {(x(k ′ ), . . . , x(k))}

Appendix 2

Therapeutic bullets found for the case study.

-

AT M golden -AT M -CHK2 golden -HRR -AT M golden -ssDN ARP A -AT M golden -BRCA1 -AT M golden -M RN -AT M golden -F AN 1 -AT M golden -ICL -DSB golden -F Acore -AT M golden -U SP 1 -AT M golden -AT M -H2AX golden -ADD -AT M golden -RAD51 -AT M golden -XP F -AT M golden -F AN CM -AT M golden -F AN CD1N -AT M golden -AT M -CHK1 golden -ICL -AT M golden -AT M -p53 golden -F AN CJBRCA1 -AT M golden -F AN CD2I -AT M golden -ICL -F AN CD1N -AT M golden -ICL -F Acore -DSB golden -BRCA1 -U SP 1 -AT M golden -BRCA1 -ssDN ARP A -AT M golden -BRCA1 -AT M -CHK1 golden -ADD -AT M -H2AX golden -F AN 1 -M RN -AT M golden -AT M -CHK2 -H2AX golden -ICL -DSB -M RN golden -XP F -M RN -AT M golden -F Acore -F AN CD2I -AT M golden -F AN CM -AT M -CHK2 golden -RAD51 -AT M -p53 golden -ICL -ssDN ARP A -AT M golden -F AN CM -AT R -AT M golden -RAD51 -AT M -H2AX golden -ADD -F AN CD1N -AT M golden -ICL -U SP 1 -AT M golden -F AN CM -M RN -AT R golden -M RN -U SP 1 -AT M golden -F AN 1 -HRR -AT M golden -BRCA1 -AT M -H2AX golden -F AN CJBRCA1 -ADD -AT M golden -M RN -ssDN ARP A -AT M golden -F Acore -ssDN ARP A -AT M golden -F Acore -F AN CD1N -AT M golden -F AN CD2I -BRCA1 -AT M golden -ADD -M RN -AT M golden -AT M -p53 -CHK2 golden -RAD51 -AT M -CHK2 golden -F AN CM -AT M -H2AX golden -ADD -P CN AT LS -AT M golden -F AN CJBRCA1 -AT M -p53 golden -F AN CM -M RN -AT M golden -F AN CJBRCA1 -AT M -CHK2 golden -F AN CD2I -U SP 1 -AT M golden -ADD -AT M -CHK2 golden -F AN CD2I -F AN CD1N -AT M golden -M RN -HRR -AT M golden -ICL -DSB -U SP 1 golden -F AN 1 -F AN CD1N -AT M golden -F AN 1 -AT M -H2AX golden -F AN CJBRCA1 -F AN 1 -AT M golden -ssDN ARP A -AT M -H2AX golden -AT M -CHK1 -CHK2 golden -ADD -HRR -AT M golden -AT M -p53 -CHK1 golden -F Acore -AT M -H2AX golden -F AN CD2I -AT M -CHK2 golden -F AN 1 -RAD51 -AT M golden -F AN CD2I -RAD51 -AT M golden -F AN CJBRCA1 -XP F -AT M golden -ICL -F AN CJBRCA1 -DSB golden -ssDN ARP A -HRR -AT M golden -M RN -BRCA1 -AT M golden -F AN CM -F AN 1 -AT M golden -ssDN ARP A -AT M -p53 golden -F AN 1 -AT M -CHK2 golden -F AN CD2I -ssDN ARP A -AT M golden -F AN CD2I -F AN 1 -AT M golden -XP F -HRR -AT M golden -F AN 1 -BRCA1 -AT M golden -ADD -AT M -CHK1 golden -F Acore -HRR -AT M golden -XP F -AT M -CHK1 golden -ADD -BRCA1 -AT M golden -ICL -F AN 1 -DSB golden -ADD -AT M -p53 golden -ICL -M U S81 -AT M golden -F Acore -RAD51 -AT M golden -AT M -CHK1 -

Appendix 3

Therapeutic bullets found for the case study using the new criterion.

Figure 1 -

 1 Figure1-Graphical representation of the example network adapted from[START_REF] Faure | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]. CDKs (cyclin-dependent kinases) are the catalytic partners of cyclins and, in this model, are not explicitly shown since the activity of CDK-cyclin complexes essentially depends on cyclins. Furthermore, the inhibition of E2F by Rb is modeled by opposing Rb to the effects of E2F on its targets. The same applies to the inhibition of CycE and CycA by p27. For a complete description of the model, see[START_REF] Faure | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]. CycD: CDK4/6-cyclin D complex, input of the model, initiates the cell cycle, activated by positive signals such as growth factors; CycE: CDK2-cyclin E complex; CycA: CDK2-cyclin A complex; CycB: CDK1-cyclin B complex; Rb: retinoblastoma protein, a tumor suppressor; E2F: a family of transcription factors divided into activator and repressor members, in this model E2F represents the activator members; p27: p27/Kip1, a CKI (CDK inhibitor); Cdc20: an APC (Anaphase Promoting Complex, an E3 ubiquitin ligase) activator; Cdh1: an APC activator; UbcH10: an E2 ubiquitin conjugating enzyme.

  a physio3 , a physio4 , a physio5 , a physio6 } where

3 k = 1 4x(k) = x 0 5 while true do 6 x 7 for targ i ∈ c targ do 8 for v j ∈ V do 9 if v j = targ i then 10 x 14 if 18 else 19 a i .f req = 1 20A

 3156789101418191 function A = compute_A(f , c targ , c moda , D, V )1 A = {} 2 for x 0 ∈ D do (k + 1) = f (x(k)) j (k + 1) = moda i ∃k ′ ∈ [[1, k]] : x(k + 1) = x(k ′ ) then 15 a i .seq = (x(k ′ ), . . . , x(k)) 16 if ∃a j ∈ A : a i .seq = a j .seq then 17 a j .f req = a j .f req + 1 = A ∪ {a i } 27 for a ∈ A do28 a.f req = a.f req • 100/card D 29 end for 30 return A end function

20 H 30 C targ = {} 31 C 33 generate c targ / ∈ C targ 34 C 35 end while 36 while card C moda < max r moda do 37 generate c moda / ∈ C moda 38 C 39 end while 40 for c targ ∈ C targ do 41 for c moda ∈ C moda do 42 H = {} 43 A patho = {} 44 for x 0 ∈ D do 45 k = 1

 203031333435363738394041424344451 = H ∪ {(x(1), . . . , x(k))} 21 end for 22 return A physio 23 prompt r min , r max , max targ , max moda 24 r max = min(r max , n) 25 golden_set = {} 26 silver_set = {} 27 for r ∈ [[r min , r max ]] do 28 max r targ = min(max targ , n!/(r! • (nr)!)) 29 max r moda = min(max moda , 2 r ) moda = {} 32 while card C targ < max r targ do targ = C targ ∪ {c targ } moda = C moda ∪ {c moda }

  www.gnu.org/software/gcc/fortran/ 2 https://github.com/ each of them attracting 50% of the x ∈ S under f physio . Then, A physio = {a 1 , a 2 } and corresponds to the results obtained by Faure et al. In terms of phenotypes, a 1 corresponds to cell cycle whereas a 2 corresponds to quiescence.

  a 1 and a 2 attracting respectively 29.5% and 70.5% of the x ∈ D under f patho2

  One of the main limitations of Boolean networks is that variables can take only two values, which can be quite simplistic. Depending on what variables model, such as activity level of enzymes or abundance of gene products, considering more than two possible levels should enable models to be more realistic. Without leaving the logic-based modeling formalism, one solution is to extend Boolean logic to multivalued logic[START_REF] Rescher | Many-valued logic[END_REF]. As with Boolean logic, variables of multivalued logic are discrete, their value belonging to [[0; 1]] where 0 means false and 1 means true. With Boolean logic, only 0 and 1 can be used to valuate variables. With multivalued logic, an arbitrary finite number h of values in [[0; 1]] can be used.

	5.1.1 Introduction

  r max ]] do test = compute_A(f patho , c targ , c moda , D, V ) 25 if A test ⊆ A physio ∪ A patho then 26 cover test = compute_cover(A physio , A test )

	20	C moda = C moda ∪ {c moda }
	21	end while
	22	for c targ ∈ C targ do
	23	for c moda ∈ C moda do
	27	if cover test > cover patho then
	28	T = T ∪ {(c targ , c moda )}
	29	end if
	30	end if
	31	end for
	32
	12	C targ = {}
	13	C moda = {}
	14	while card C targ < min(max targ , n!/(r! • (n -r)!)) do
	15	generate c targ / ∈ C targ
	16	C targ = C targ ∪ {c targ }
	17	end while
	18	while card C moda < min(max moda , h r ) do
	19	generate c moda / ∈ C moda

24

A

x(k) = x 0 while true do if ∃w ∈ H : x(k) ∈ w then break end if

end for end for end for return golden_set, silver_set