BACKGROUND
the emerging field of systems biology aims to gain high-level understanding of complex biological systems by studying the structure and dynamics of cellular functions  <cit> .
the construction and analysis of models describing biological processes are at the core of this new approach.
computer science provides methods for simulating and analyzing these models.
this allows performing in silico experiments that lead to new predictions and thus complement traditional experimental approaches.
considering biology from a systems point of view also allows the use of methods designed for the engineering of artificial systems.
a wide array of formalisms originating from computer science have already been used to model biological systems.
these include process calculi  <cit> , statecharts and live sequence charts  <cit> , variants of petri nets  <cit> , and hybrid models  <cit> .
recently, the term executable biology has been proposed for models that consist of an algorithm which mimics the behavior of a studied biological system  <cit> .
executable models of biological systems can be analyzed using formal verification methods, which were originally designed for the validation of computer systems.
two possible approaches are testing and model checking.
model checking  <cit>  consists of an exhaustive exploration of the state-space of a given system in order to verify that it always adheres to a set of requirements.
in the case of biological systems, requirements are derived from experimental observations.
the consistency of a putative executable model can thus be tested with respect to the data on which the model is based  <cit>  .
in order to construct executable models of biological pathways, it is necessary to precisely define the interactions between biochemical components.
one of the widely used formalism for representing interactions in biology are boolean networks  <cit> .
a boolean network consists of a set of boolean variables representing the activity of biochemical components such as proteins or genes.
a boolean function is associated to every component.
this function updates the value of the corresponding component at every step of the execution of a boolean network.
in the context of biological pathways, a boolean network can be built from a graph representing biochemical interactions.
every node of this graph represents a biochemical component.
the arcs of the graph are used to represent the interactions between components.
their weight can be positive  or negative .
this graph is used to define boolean functions for every component.
the function sets the value of a component to true if the amount of activation is higher than the amount of inhibition.
else the value of the component is set to false.
boolean networks therefore provide a highly abstracted view of protein activity or gene expression.
this level of abstraction often corresponds to the available qualitative experimental data.
furthermore, the interactions between components are represented in a very similar way in the diagrammatic models that are widely used by biologists.
boolean networks therefore appear to be well suited for the qualitative analysis of biological processes.
they have been used recently for modeling a wide array of biological systems, such as the arabidopsis thaliana flower morphogenesis  <cit> , the influence of physical stimuli on cell fate  <cit> , the segment polarity gene network in drosophila melanogaster  <cit>  and the yeast cell cycle  <cit> .
boolean networks have been extended in several ways, including probabilistic boolean networks  <cit>  and the so-called generalized logical approach using multivalued levels of concentration  <cit> .
the generalized logical approach has been combined with formal verification methods  <cit> , and constraint programming was used for finding steady states of large networks  <cit> .
in this work, we propose an extension of boolean networks called qualitative networks.
we use this framework to construct an executable model which is analyzed using formal verification methods.
we introduce two changes to the boolean networks framework.
first, we use discrete variables instead of boolean variables for representing gene expression and protein activity levels.
second, we capture a richer set of interactions between components than only activation and inhibition.
any function of the current state of the system can be used to define the target towards which a component will move at the next step.
at each step, the level of a component changes towards this target by at most one level.
we use simulation and model checking to verify that the implementation of the qualitative network adheres to a set of specifications derived from the laboratory experimental observations.
if necessary, we can incorporate new hypotheses into the model .
once these are tested in-silico, they should also be validated experimentally.
we propose that the stable states of the biological system correspond to the infinitely visited states of the executable model.
we introduce an efficient symbolic algorithm for verifying if all infinitely visited states of a qualitative network satisfy the specifications.
we consider the performance of this algorithm for models of different size and complexity, and show that it can be used to analyze large and complex qualitative networks.
we illustrate the ability of qualitative networks to provide a valid qualitative representation of a studied biological system using a model of the incoherent type- <dig> feed-forward loop in the e. coli gal system.
our approach to the modeling of signaling pathways allows the detection of gaps in the mechanistic understanding of the studied process, and suggests new biochemical interactions at a similar level of abstraction to the one observed in diagrammatic models.
we illustrate this approach by a multicellular model of the interaction between the notch and the wnt pathways in mammalian skin cells.
RESULTS
qualitative networks
in this study we propose an approach for constructing executable models of signaling pathways which is based on qualitative networks, an extension of boolean networks.
in this section, we first explain how biological pathways can be modeled using boolean networks and then define qualitative networks.
interactions in biological pathways can be represented as an edge-weighted interaction graph g .
a node vi ∈ v represents a biological component such as a gene or a protein.
an edge eij ∈ e from node vi to node vj has a weight αij, corresponding to the effect of the component represented by vi on the component represented by vj.
activation is specified by a positive value of αij, inhibition by a negative value.
this graph therefore represents the type of interaction between components as well as their strength in a similar way to the diagrammatic models used by biologists.
li et al.
propose the use of boolean networks for studying the interactions defined in such a graph  <cit> .
a boolean network b  consists of a set of components c and a list of boolean functions f. each component ci ∈ c has a state, which can take two values: ci =  <dig> if the component is active and ci =  <dig> if the component is absent or inactive.
we use k to represent the number of components .
the state s of a boolean network corresponds to the vector  ∈ { <dig>  1}k of the states of all components of the boolean network.
a boolean function fi ∈ f is assigned to every component.
this function maps the current state of the system to a boolean value representing the next state of the corresponding component .
time in boolean networks is represented by a discrete variable t. we use the notations ci  and s  to represent the state of a component, respectively the state of the whole system, at time t. the value of each protein at the next time step is computed by applying the corresponding boolean function to the current state of the system:
ci  = fi )
given an interaction graph g , a boolean network b  can be obtained as follows: each node of v corresponds to a component in c. the boolean function should reproduce the following behavior for every component: if the cumulative activation is higher than the cumulative inhibition, the state of the component is set to active  and, vice-versa, if the inhibition is higher than the activation, it is set to inactive .
hence, the boolean function corresponding to each component is defined as follows:
fi={0if∑cj∈sαjicj<01if∑cj∈sαjicj>0ciif∑cj∈sαjicj=0
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgmbgzdawgaawcbagaemyaakgabeaakiabcicaoiabdohazjabcmcapiabg2da9maaceqabaqbaeqabmgaaaqaaiabicdawaqaaiabbmgapjabbagamnaaqafabaaccigae8xsde2aasbaasqaaiabdqgaqjabdmgapbqabagccqwgjbwydawgaawcbagaemoaaogabeaaaeaacqwgjbwydawgaaadbagaemoaaogabeaaliabgigiolabdohazbqab0gaeyyeiuoakiabgyda8iabicdawaqaaiabigdaxaqaaiabbmgapjabbagamnaaqafabagae8xsde2aasbaasqaaiabdqgaqjabdmgapbqabagccqwgjbwydawgaawcbagaemoaaogabeaaaeaacqwgjbwydawgaaadbagaemoaaogabeaaliabgigiolabdohazbqab0gaeyyeiuoakiabg6da+iabicdawaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaagcbagaeeyaakmaeeozay2aaabuaeaacqwfxoqydawgaawcbagaemoaaomaemyaakgabeaakiabdogajnaabaaaleaacqwgqbgaaeqaaaqaaiabdogajnaabaaameaacqwgqbgaaeqaasgaeyici4saem4camhabeqdcqghris5aogaeyypa0jaegimaadaaagaay5eaaaaaa@759c@
we extend the boolean network framework by using discrete domains to represent the state of a component and by allowing a rich set of interaction between components.
the state of an active component can take several values, rather than being simply on.
these values correspond to qualitative levels of gene expression or protein activity, for example, low, medium, high.
this representation therefore is a middle ground between the boolean abstraction and quantitative models.
the use of discrete domains requires to modify the way the updates of the network are defined.
we introduce a target function that computes the discrete level towards which a component moves in the next step.
the target function is not restricted to modeling activation and inhibition.
every function of the current state of the system can be used as target function.
this makes it possible to model interactions that cannot be represented in interaction graphs.
one such example is complexation, where the amount of produced protein is bounded by the availability of the proteins required for the interaction.
this situation can be modeled by using the minimum of the level of these proteins as a target function.
more elaborate target functions can also be used to represent situations in which the effect of a specific activating or inhibiting edge is dependent on the level of additional proteins.
we call this extended framework qualitative network.
in the following paragraphs, we give a detailed description of qualitative networks.
a qualitative network q  consists of a set of components c and a list of target functions t. the state of a component ci ∈ c is an integer value between  <dig> and n. these values represent qualitative levels of the components with  <dig> the lowest and n the highest possible level for that component.
for each component ci, there exists a target function targeti ∈ t which specifies the level towards which the value of the component should move.
this function is computed on the state of the whole system .
at each step, the value of the state of a component can only increase or decrease by a single level.
the value of each component at the next time step is therefore computed as follows:
ci={ci−1if targeti)<cici+1if targeti)>ciciif targeti)=ci
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabgucariabigdaxiabcmcapiabg2da9maaceqabaqbaeaabmgaaaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaogaeiikagiaemidaqnaeiykakiaeyoei0iaegymaedabagaeeyaakmaeeozaymaeeiiaaiaemidaqhcbigae8xyaemae8ncainae83zacmaemyzaumaemidaq3aasbaasqaaiabdmgapbqabagccqggoaakcqwgtbwucqggoaakcqwg0badcqggpaqkcqggpaqkcqgh8aapcqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabcmcapaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaogaeiikagiaemidaqnaeiykakiaey4kasiaegymaedabagaeeyaakmaeeozaymaeeiiaaiaemidaqnae8xyaemae8ncainae83zacmaemyzaumaemidaq3aasbaasqaaiabdmgapbqabagccqggoaakcqwgtbwucqggoaakcqwg0badcqggpaqkcqggpaqkcqgh+agpcqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabcmcapaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaogaeiikagiaemidaqnaeiykakcabagaeeyaakmaeeozaymaeeiiaaiaemidaqnae8xyaemae8ncainae83zacmaemyzaumaemidaq3aasbaasqaaiabdmgapbqabagccqggoaakcqwgtbwucqggoaakcqwg0badcqggpaqkcqggpaqkcqgh9aqpcqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabcmcapaaaaiaawuhaaiaaxmaacawljawaaewaaeaacqaiyagmaiaawicacaglpaaaaaa@9f11@
when n =  <dig>  this update function is equivalent to the update function of boolean networks .
it follows that qualitative networks are a generalization of boolean networks: the qualitative network q  is equivalent to the boolean network b .
a particular kind of target function is directly obtained from the qualitative graph describing the studied biological system.
it is used to model activation and inhibition in a similar way to the one describe for boolean networks.
in boolean networks, the sum of all contributions was sufficient to decide on the next value of a component.
here we need to obtain a target function that ranges over { <dig>  ..., n}, and for which intermediate values are meaningful.
in order to obtain the target function for component ci, we compute separately the amount of activation acti and the amount of inhibition inhi sensed by this component.
both of these values are scaled so that they are proportional to the weighted amount of activation, respectively inhibition.
when all activating components are at their highest value, then acti = n. similarly, when all inhibiting components are at their highest value then inhi = n.
acti=∑αji>0αjicj∑αji>0αjiinhi=∑αji<0αjicj∑αji<0αji
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakqaabeqaaiabdggahjabdogajjabdsha0naabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakiaeyypa0zaasaaaeaadaaeqbqaaggaciab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaogaem4yam2aasbaasqaaiabdqgaqbqabaaabagae8xsde2aasbaawqaaiabdqgaqjabdmgapbqabawccqgh+agpcqaiwaamaeqaniabgghildaakeaadaaeqbqaaiab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaaqaaiab=f7ahnaabaaameaacqwgqbgacqwgpbqaaeqaasgaeyopa4jaegimaadabeqdcqghris5aaaaaoqaaiabdmgapjabd6gaujabdigaonaabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakiaeyypa0zaasaaaeaadaaeqbqaaiab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaogaem4yam2aasbaasqaaiabdqgaqbqabaaabagae8xsde2aasbaawqaaiabdqgaqjabdmgapbqabawccqgh8aapcqaiwaamaeqaniabgghildaakeaadaaeqbqaaiab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaaqaaiab=f7ahnaabaaameaacqwgqbgacqwgpbqaaeqaasgaeyipawjaegimaadabeqdcqghris5aaaaaaaa@7aa5@
the target function is the difference between the amount of activation and the amount of inhibition.
when there is no inhibition, the output of the target function corresponds to the amount of activation.
the level of the component will not exceed the amount of the activation.
when the amount of inhibition is higher than the amount of activation, the value of the target function is zero.
we need to consider a special case when all incoming edges of a component are inhibitions.
in this case, we must assume that the mechanism leading to the presence of this component is not modeled, and that, in the absence of inhibition, the target function of this variable will return the maximal value.
we obtain the following definition of the target function:
targeti={max⁡−inhi)if max>0n−inhiif max≤0
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwg0badiegacqwfhbqycqwfybgccqwfnbwzcqwglbqzcqwg0baddawgaawcbagaemyaakgabeaakiabcicaoiabdohazjabcmcapiabg2da9maaceqabaqbaeaabigaaaqaaigbc2gatjabcggahjabciha4jabcicaoiabicdawiabcycasiabdggahjabdogajjabdsha0naabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakiaeyoei0iaemyaakmaemoba4maemiaag2aasbaasqaaiabdmgapbqabagccqggoaakcqwgzbwccqggpaqkcqggpaqkaeaacqqgpbqacqqgmbgzcqqggaaicqwftbqbcqwfhbqycqwf4baecqggoaakiigacqgfxoqydawgaawcbagaemoaaomaemyaakgabeaakiabcmcapiabg6da+iabicdawaqaaiabd6eaojabgkhitiabdmgapjabd6gaujabdigaonaabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakcabagaeeyaakmaeeozaymaeeiiaaiae8xba0mae8xyaemae8heagnaeiikagiae4xsde2aasbaasqaaiabdqgaqjabdmgapbqabagccqggpaqkcqghkjyocqaiwaamaaaacagl7baacawljagaaczcamaabmaabagaeg4mamdacagloagaayzkaaaaaa@84bd@
we use qualitative networks to define executable models of biological systems.
the list of target functions specifies how the different molecules in the system influence each other.
when creating a model based on the mechanistic understanding of a biological process, we use this information to define a set of target functions.
where necessary, we separate the representation of protein activity from the expression level of the corresponding gene.
this is useful for modeling cases in which there are interactions at the transcription level as well as protein-protein interactions.
during the analysis process, we may modify the definition of the target functions.
while qualitative networks allow for any kind of target function, we use a subset of functions that represent specific biological interactions, such as the activation/inhibition function described above.
this ensures that new hypothesized interactions between components are defined at a level where they have a clear biological meaning.
the implementation of models defined as qualitative networks is described in the methods section.
iterative model improvement
a model in general, and an executable model in particular, is meant to represent how the system that we study actually works.
thus, it is natural to expect that the model would be consistent with the data obtained from laboratory experimentation with this system.
we use simulation and model checking  <cit> , which allows to formally check all the different executions of the model against a formal specification.
by exploring all the possible states and transitions of the model we can determine whether some requirement holds for the model.
in the case that the property does not hold, the model-checking algorithm supplies a counter example: an execution of the system that does not satisfy a given requirement.
we use the information provided by this counter example to suggest modifications of the qualitative network representing the studied system.
our methodology in using qualitative networks is the iterative improvement process .
we stress that the usage of qualitative networks is not restricted to this methodology.
the advantages of multiple values, flexible target functions, and the analysis algorithms suggested below, can be used regardless of the modeling methodology chosen by the user.
the iterative improvement process consists of the following.
we first derive a set of specifications from prior laboratory experimental results.
we build an initial putative model based on the current mechanistic understanding of the studied biological process and the working hypotheses we would like to evaluate.
we then use formal verification methods to verify if the proposed model conforms with the specification.
if the model fails this test, we deduce that the model does not conform to the data.
we use counter examples, to suggest a revised model by making new hypotheses.
the new model is then evaluated again, and these steps are repeated until no contradiction is found between the model and the specifications.
at this point, we might have had to modify our hypotheses, both in discarding unverifiable assumptions and in adding assumptions that are necessary for the model to comply with the specifications.
the resulting model, and the corresponding new set of hypotheses, are then known to be consistent with the experimental results.
at this stage, it is possible to query the model to gain additional knowledge about the biological system.
this can be done either by modifying the model to evaluate or by adding specifications that could provide additional information on the behavior of the system.
hypotheses and changes in the model leading to interesting behaviors are candidates for experimental validation.
the result of such an experiment, whether confirming the information gained from the computational analysis or not, can then be added to the specification.
this creates a second, outer, improvement loop which combines computational and experimental validation.
the specifications derived from laboratory experimental results can be separated into two categories.
low-level specifications involve changes in the level of protein activity or gene expression .
high-level specifications consider changes at the cellular level, such as cell fate or susceptibility to carcinogenesis.
since we can only observe the level of individual proteins, we need to define how these relate to changes at the cellular level.
experimental results are often related to manipulation, such as gene over-expression or knockout.
these manipulation lead to changes in the behavior of the model.
the methods section provides additional details on how specifications are defined and verified.
this iterative improvement approach can be applied to highly non-deterministic models, hence the use of model checking.
qualitative networks, like boolean networks, have a deterministic behavior, and it is therefore sufficient to perform one execution per initial state, and verify if each of these executions satisfies the specifications.
if some executions do satisfy the specifications, but others do not, then the model is not constrained enough, and we need to formulate hypotheses that allow avoiding the executions which do not satisfy the specifications.
ideally, the new hypotheses would not contradict the current mechanistic understanding of the system, but rather suggest additional putative interactions between components.
if there is a specification which is not satisfied by any execution of the model, then there is a contradiction between the mechanistic understanding of the system and the experimental data.
in this case, the modifications needed for the model to satisfy the requirement will suggest alternative mechanistic models of the studied system that are more consistent with the experimental data.
in the following sections, we describe two additional methods for analyzing the model during the improvement process: the use of non-determinism on parts of the model, and the efficient computation of all infinitely visited states.
use of non-determinism
in the case where every interaction in the model is precisely defined, enumerative networks are deterministic.
non-determinism can, however, be used to represent unknown interactions.
in particular, rather than studying a complete multicellular model, the model of a single cell can be studied by representing inter-cellular interactions non-deterministically.
the level of proteins that are external to the cell are modeled by non-deterministic variables.
these variables cannot jump arbitrarily from one value to another value, but only change by one level per time step.
at each time step, the variable can either increase, decrease or remain constant.
this is sufficient to capture any possible behavior since updates of components are bounded by equation  <dig>  the set of possible behaviors of the state of this protein is therefore a superset of the possible behaviors for any possible target function.
instead of simulation, we use model checking, which explores all possible executions of a non-deterministic system.
this leads to the analysis of a superset of the possible executions of any model with more constrained inter-cellular interactions.
if no execution of the system satisfies the specification, then it is impossible that any model including this cell can satisfy the specifications, no matter how the inter-cellular interactions are defined.
if some executions satisfy the specifications, but other do not, then it is not possible to predict the behavior of the cell with more constrained inter-cellular interactions.
this approach is therefore useful for finding out if the improvement should be performed at the level of the inter-cellular communication or at the intracellular level.
finding infinitely visited states
we propose an efficient symbolic algorithm for computing the infinitely visited states of a qualitative network.
a state s is said to be infinitely visited if there is an execution of the network in which the state appears infinitely often.
in a biological context, the infinitely visited states correspond to the stable states of the system.
in the context of specification-based analysis of the model, we are therefore interested in knowing if the specification holds on all infinitely visited states.
this captures the fact that, as opposed to computerized systems, biological system do not have an initial state, but rather are continuously behaving.
after cell division, for example, the two resulting cells would have proteins levels that are at a qualitative level similar to the levels in the original cell.
states that are not infinitely visited can thus be considered as unstable states.
if the model is in such a state, it will evolve towards an infinitely visited state in a bounded number of steps, and then remain in infinitely visited states for an indefinite time.
therefore, we do not insist that the specification holds on states that are not infinitely visited.
analysis of boolean networks also include studies of the attractors of the model.
in discrete, deterministic models, attractors are loops of one or more states that are visited infinitely often.
questions of interest include: does the system end up in one or several attractors, and what initial conditions lead to which attractor.
this analysis requires exploring the executions starting from all initial states.
the number of possible initial states is exponential in the number of variables, and therefore we try to prune them by using the available information on the system.
in case that no biologically meaningful definition of the set of initial state exists, it is necessary to consider all states as being initial.
in this case, enumerative exploration is practically intractable even for relatively small models.
in the methods section, we describe a symbolic algorithm for finding all the attractors of a model, based on compact representation of sets of states using binary decision diagrams   <cit> .
the algorithm we propose uses the structure of multicellular qualitative networks by interleaving composition and computation of infinitely visited states.
we first consider a partial model consisting only of the first cell.
the behavior of components in neighbor cells are abstracted by allowing them to change non-deterministically.
we compute the set of infinitely visited states of this partial model.
this set contains the projection of any infinitely visited state of the complete model onto the variables of the partial model.
we use this information to compute the infinitely visited states of a partial model consisting of the first two cells.
this process is repeated to obtain the complete model and the set of infinitely visited states.
we further optimize this algorithm by using partitioned transition relations  <cit> .
based on the target functions list of the qualitative network, we define a transition relation for a single cell.
we then use it repeatedly to compute the next state of the multicellular system, rather than using a single transition relation for the whole model.
the methods section contains a detailed description of this algorithm.
modeling network motifs using qualitative networks
transcription regulation networks contain patterns that occur significantly more often in a real network than in a random network with the same characteristics.
these patterns are called network motifs  <cit> .
the multiple values used by qualitative networks may be required to accurately model the qualitative behavior of common network motifs in a way that would be difficult, if not impossible, to achieve with boolean values.
we illustrate this ability by considering the case of the incoherent type- <dig> feed-forward loop  in e. coli.
in this motif, transcription factor x activates transcription factor y and gene z, and transcription factor y represses z. the response time is defined as the time needed for z to increase halfway to its steady-state level.
compared to simple regulation , the i1-ffl has a shorter response time to the activation of x. when x is activated, the level of both y and z start to increase.
once the level of y exceeds the repression threshold for the z promoter, it will start repressing z. as a consequence, the level of z starts decreasing towards a low steady-state level.
the response time is thus shorter, because the steady-state level of z is lower than for simple activation, but its initial transcription rate is the same in both cases.
this behavior has been observed both using quantitative models  <cit>  and experimentally in the gal system in e. coli  <cit> .
in the gal system , an i1-ffl links the activator crp to the galactose utilization operon galetk.
when there is glucose starvation, crp gets activated.
crp then activates both galetk and gals, and gals inhibits galetk, hence the i1-ffl motif.
in addition, β-d-galactose or d-fucose  unbind gals from the galetk promoter region, thus cancelling its inhibitory effect.
finally, gals negatively auto-regulates its own expression.
we construct a qualitative network model of this system based on these interactions.
this model needs to satisfy two requirements to be realistic.
first, the activation of crp must lead to the activation of galetk independently of the level of galsinh.
second, the time needed by galetk to increase halfway to its steady-state needs to be shorter when galsinh is absent than when it is present.
figure 2b shows the evolution over time of the level of galetk after crp is activated for both case.
this result is qualitatively similar to the experimental measurements of the galetk expression performed by mangan et al.
<cit> .
qualitative networks are able to provide realistic qualitative approximations of the behavior of further common network motifs, such as other types of feed-forward loops and regulated-feedback motifs that appear in developmental transcription networks  <cit> .
case study: a model of the interaction between the notch and wnt signaling pathways
we apply our modeling approach to a model representing the crosstalk between the notch and wnt pathways in mammalian keratinocytes.
both pathways play a key role in the control of cell proliferation and differentiation and have been linked to several types of cancer.
we first present the biological background, then introduce the qualitative network model for a single cell and for a model with five cell representing the multiple layers of the mammalian epidermis.
we define a set of specifications derived from laboratory experimental perturbations on these pathways, and use them to analyze the possible executions of the model.
crosstalk between the notch and wnt signaling pathways
notch signaling plays a key role in the control of the development of various tissues, with the particularity that, depending on the context, it can either induce differentiation or maintain cells in a proliferating state .
the notch receptor protein is located across the cell membrane and is composed of two parts – an extracellular component and a cytoplasmic part.
the extracellular component of notch can bind to a ligand expressed on the membrane of a neighboring cell.
in mammals, there are four types of notch proteins, three types of delta-like ligands and two types of serrate-like ligands, which are called jagged.
the binding of a ligand to the notch receptor leads to a proteolytic cleavage of the two parts of the notch protein, thus releasing its cytoplasmic part – the notch intracellular domain .
notch-ic then travels to the nucleus, where it binds to the transcription factor csl.
in the absence of notch-ic, csl acts as a transcription inhibitor, but the presence of notch-ic converts it to a transcription activator, thus leading to the expression of its transcriptional targets.
the wnt signaling pathway has a critical regulating role in stem cells, and has been associated with cancer in several tissues .
the central player in the canonical wnt signaling pathway  is the protein β-catenin.
in the absence of wnt signaling, the newly formed β-catenin is found in a degradation complex composed of the proteins adenomous polyposis coli , the scaffolding protein axin and the glycogen-synthase kinase-3β .
the gsk3β phosphorilates β-catenin, which is then degradated by proteasomes.
wnt is a short range signaling protein which interacts with a serpentine receptor of the frizzled family.
the consequence of this interaction are not completely understood, but it appears that an axin-binding molecule, dishevelled , gets phosphorylated and starts inhibiting the degradation of β-catenin.
this β-catenin can therefore accumulate in the cytoplasm and travel to the nucleus, where it interacts with the tcf/lef dna-binding proteins.
in the absence of β-catenin, a co-repressor like groucho binds to tcf/lef, thus blocking the transcription of its downstream target genes.
when interacting with β-catenin, tcf/lef becomes a transcription activator, thus leading to the activation of the downstream target genes of the wnt signaling pathway.
the outermost layer of the human skin, the epidermis, can be divided into several layers: the basal layer, which is the closest to the underlying dermis, the suprabasal layer, and the cornified layer, which is mainly composed of dead cells .
in normal tissue, most proliferating keratinocytes are located in the basal layer.
they migrate to the suprabasal layer when they have initiated terminal differentiation and then stop proliferating.
notch is expressed in the human epidermis, at a higher level in the suprabasal layers than in the basal layer.
the jagged ligand is co-expressed with notch in the skin and has been proposed as a positive feedback mechanism for sustained notch signaling triggering differentiation  <cit> .
the interaction between the notch and the wnt pathways are involved in the transition from a proliferating to a differentiated state in keratinocytes.
in vitro experiments have shown that notch signaling induces terminal differentiation in both murine and human skin  <cit> .
the tissue specific ablation of the notchl receptor in mouse epidermis leads to hyperplasia in the epidermis, the development of tumors and an increased susceptibility to chemical carcinogens – thus indicating that notch acts as a tumor suppressor in mammal skin  <cit> .
notch was suggested to interact with the wnt pathway through its direct transcriptional target p <dig>  which is a negative transcription regulator of wnt  <cit> .
model construction
we build a model of the notch and wnt pathways in the epidermal layer of mammalian skin.
our model is composed of five identical keratinocytes, each of them representing one layer of the epidermis.
protein activity and gene expression are represented by variables that can take four values: off, low, medium and high.
all activation and inhibition interactions have an equal weight .
figure  <dig> represents the two pathways in a single cell.
except for the intracellular part of notch , the target function of all other components is the difference between activation and inhibition, as defined by equation  <dig>  the target function of the notch-ic component is the minimum between the level of notch receptor on the cell membrane and the level of ligand on the membrane of the neighboring cells.
the complete model consists of a single row of five cells, which are named cell <dig> to cell <dig>  neighboring cells are connected both through the wnt and the notch pathways .
cell <dig>  the leftmost cell on the visualization of the model, represents the lowest layer of keratinocytes of the epidermis, which means that the dermis is to its immediate left.
this cell should adopt a proliferating fate.
the upper skin, to which the keratinocytes migrate before they die, is represented on the right side of the model.
we consider that the two immediate neighbors contribute to the level of ligand sensed by a cell.
the level of ligand sensed by a cell  is activated by the level of ligand on the membrane of its immediate neighbor, and is thus maximal when the level of ligand in both neighbors is high.
the level of jagged on both extremities of the model is constant at medium.
similarly, we connect the level of wnt produced by a cell to the level of wnt sensed by the nearby cells.
we consider the case in which a cell does not sense wnt emitted by itself.
in this case we can use the same double activation scheme as for the ligands.
on the left side of the model, the level of wnt is fixed to high since the dermis is known to emit wnt signaling proteins.
in contrast, the level on the right side is fixed to low since wnt signaling does not occur in the upper layers of the skin.
specifications
in order to be able to define high level specifications, we need to map changes in the level of protein activity to a certain cell fate.
we consider the balance between the downstream target genes of wnt signaling  and the downstream target genes of notch signaling .
the downstream targets of the wnt signaling are known to maintain the cells in a proliferating state, whereas the downstream targets of notch initiate differentiation.
therefore, we make the assumption that a cell in which gt <dig> > gt <dig> is proliferating, whereas a cell in which gt <dig> < gt <dig> is differentiated.
furthermore, since an increase in the level of gt <dig> or a decrease in the level of gt <dig> makes a cell more likely to be proliferating, we assume that this change also makes the cell more sensitive to chemically-induced carcinogenesis.
for experiments involving the knockout or the over-expression of a gene, we compare the level of components when the wild type model is in a stable state to the level of the proteins when the altered model reaches a stable state again.
the first high level specification  defines the wild type state, and can be split into two parts: cell <dig> is proliferating , and cell4– <dig> are differentiated .
h <dig> : gt <dig> > gt <dig> in cell1
h <dig> : gt <dig> < gt <dig> in cell4–5
we do not use cell2– <dig> in these specification, but if a cell has adopted a differentiated fate, cells in a higher sub-layer  cannot be proliferating.
we formulate this as an additional requirement:
h <dig>  gt <dig> = gt <dig> in celli ⇒ ∀j > i, gt <dig> ≤ gt <dig> in cellj.
h <dig>  gt <dig> < gt <dig> in celli ⇒ ∀j > i, gt <dig> < gt <dig> in cellj.
we derive further specifications from the experiments performed by nicolas et al.
on mouse keratinocytes  <cit> .
they note that notch knockout leads to increased proliferation, which we translate into a requirement that cell <dig> is also proliferating.
h3: notch = ko ⇒ gt <dig> > gt <dig> in cell4
after several months, mice in which notch is knocked out begin to develop basal cell-like carcinoma.
they are also significantly more sensitive to chemically induced carcinogenesis.
in order to relate these long-term changes to variations at the protein level, we assume that the short term effect of the notch knockout is an increase of gt <dig> in several cells, or a decrease of gt <dig> in several cells, or both.
h4: notch = ko ⇒ gt <dig> increases in cell1– <dig> or gt <dig> decreases in cell1–5
the impact of p <dig> knockout was also considered.
mice with p <dig> deficiency were more sensitive to chemically induced carcinogenesis, but did not develop tumors on their own.
we relate this result to the level of gt <dig> and gt <dig> in a similar way to notch knockout.
h5: p <dig> = ko ⇒ gt <dig> increases in cell1– <dig> or gt <dig> decreases in cell1–5
the work by nicolas et.al.
on mice keratinocytes  <cit> , as well as the work of devgan et.al.
<cit>  on the link between notch and wnt through p <dig>  also leads to several low level specifications.
these specifications were directly considered when constructing the model and are not detailed here.
analysis of the model
we compute the set of infinitely visited states of the model described above.
we obtain a total of  <dig> states.
all of these states adhere to the specifications h <dig> and h <dig>  these states, however, differ as far as the fate of cell <dig> is concerned.
this cell is either proliferating , differentiated  or in transition from proliferating to differentiated .
several of these attractors consist of more than one state, and the relation between gt <dig> and gt <dig> in cell <dig> can change between the different states of a single attractor.
we can therefore conclude that, in this model, cell <dig> is no longer receiving constant signals for proliferation, but has not yet committed to terminal differentiation.
after the model reaches a stable state, modifying it to knockout notch or p <dig> allows to verify that the model also satisfies specifications h <dig>  respectively h <dig>
further insights can be gained by studying variations of the model.
we consider the hypothesis that notch-ic activates the transcription of β-catenin, thus partially compensating its inhibitory influence on β-catenin through the wnt pathway.
this hypothesis can be rejected by analyzing a single cell model with non-deterministic levels for both external wnt and sensed ligand .
no execution of this model can satisfy requirement h <dig>   for all the states of the execution that occur infinitely often.
this means that in any multicellular model with this interaction, no cell could be proliferating.
we can therefore conclude that some other mechanism activating β-catenin must exist.
notch-ic may play a role in the activation of β-catenin, but cannot be the only activating protein.
we also consider the hypothesis that notch signaling uses the delta ligand instead of the jagged ligand.
while there is a positive feedback mechanism from notch-ic to jagged, notch-ic inhibits the expression of the delta ligand on the cell surface.
we find that such a model cannot reproduce the wild type behavior .
by analyzing the counterexample traces, we notice that due to the inhibition of the ligand, sustained notch signaling cannot be achieved in the higher layers of the epidermis.
this is consistent with the experimental analysis of the level of delta in the epidermis showing that delta is expressed in the basal layer only  <cit> .
in their analysis, lowell et al.
describe stem-cell clusters with local concentration of delta.
in order to model and analyze local effect such as cell clusters, it would be necessary to consider significantly larger multidimensional models.
performance analysis
in this section, we describe the performance of the symbolic algorithm used to compute all infinitely visited states of a qualitative network.
we first consider variations of the notch/wnt model described above.
we then use an arbitrary model to show that the method can be applied to models of more complex interactions between components.
finally, we show that our algorithm performs better than simpler symbolic algorithms.
the methods section gives additional details about the settings used for these analyses.
the computation of all infinitely visited states of the five-cells notch/wnt model takes  <dig> minutes.
this model has  <dig> ≈  <dig> initial states.
this is a significantly larger number than previously studied by boolean networks .
in order to evaluate the performance of our symbolic algorithm on larger models, we extend the notch/wnt model up to  <dig> cells.
we add cells both to the basal layer  and to the suprabasal layer .
as in the five-cells model, the level of notch receptor is set to off in the basal layer, and to high in the suprabasal layer.
while the original model has only one cell with an intermediate level of notch receptor, we introduce a second cell in the larger model in order to obtain a gradient.
these models have been developed for the purpose of performance analysis only, and do not provide additional biological insights.
table  <dig> shows the time needed to compute the set of all infinitely visited states for models with  <dig> to  <dig> cells.
the results illustrate that we can compute the set of infinitely visited states of a model which has more than  <dig> initial states.
it can be observed that adding a cell in the basal layer has a significantly smaller impact on the performance than adding a cell in the suprabasal layer.
this is due to the fact that the possible behaviors of the cells in the basal layer are very limited due to the lack of notch signaling.
in general, the performance of symbolic algorithms is dependent on how compact the representation of the sets of states and the transition relation are.
this complexity is directly dependent on the structure of the model.
in the case of qualitative networks, choosing a different combination of target functions has an impact on the performance even though the number of components stays the same.
we illustrate this by constructing an arbitrary qualitative network , which has no biological meaning except that we only use the activation/inhibition target function.
this model contains the same number of components per cell as the notch/wnt model.
the interactions between the components are more complex.
in order to compare the complexity of the two models, we consider the size of the bdd used to represent the transition relation of a single cell.
this size, measured as the number of nodes in the bdd, relates more or less to the required memory.
the arbitrary qualitative network is represented using 279' <dig> bdd nodes, compared to 39' <dig> bdd nodes for the notch/wnt model.
this increase in complexity is reflected in the performance of the infinitely visited states computation.
when considering only three cells this computation takes  <dig> seconds for the notch/wnt model, and  <dig> minutes for the arbitrary qualitative network.
our algorithm interleaves the composition of a multi-cellular model with the computation of infinitely visited states.
this approach outperforms symbolic algorithms that directly compute the set of infinitely visited states of the complete model.
the application of the direct algorithm to models of the notch/wnt pathway with four or more cells fails, because the symbolic representations of the sets of states becomes too large.
on the three-cell model, our algorithm needs  <dig> seconds and the direct algorithm  <dig> seconds.
the advantage of first computing the infinitely visited states on partial models is that the sets of potentially infinitely visited states is smaller than in the direct algorithm.
this is illustrated in table  <dig>  which provides details about the application of our method to the ten-cells notch/wnt model.
for each partial model, we indicate the initial set of states we start from, as well as the set of infinitely visited states of the partial model after the computation.
in particular, it can be noted that the computation of the infinitely visited states of the complete model starts from a set of approximately  <dig> states.
an algorithm that directly computes all infinitely visited states would start from the set of all states of the system .
the symbolic approach we propose is able to compute the infinitely visited states for large multicellular qualitative networks and to handle complex models.
by tailoring the algorithm specifically for qualitative networks, we can analyze larger models than by directly computing the infinitely visited states.
discussion
in this work, we extend boolean networks to qualitative networks; networks that allow variables to take discrete values representing qualitative levels of protein activity or gene expression.
by introducing a target function, we allow a rich set of interactions between components.
we propose a method for analyzing qualitative networks, that uses formal verification of specifications derived from experimental observations for iteratively improving the model.
in order to efficiently analyze the properties of infinitely visited states of large models, we introduce a symbolic algorithm which builds on the structure of our modeling method.
we apply this method to build and analyze a multicellular model of the interactions between the wnt and notch pathways in mammalian keratinocytes.
we show that the proposed symbolic algorithm can be applied to large and complex multicellular models.
model analysis approach
our modeling approach is aimed at using formal verification methods for the analysis of signaling networks in an iterative improvement process which is based on specifications that are directly derived from experimental results.
this work therefore combines formal verification methods with the widely used formalism of boolean network, and its extension to qualitative domains.
since boolean networks are a particular case of qualitative networks, the iterative improvement approach can be directly applied to boolean networks as well.
our method allows to build a model at a level of abstraction that is similar to the one commonly used in cell biology, and to analyze it using specifications that are directly derived from experimental studies.
this approach is therefore useful for biologists to identify gaps between the hypothesized mechanistic description of a pathway of interest, and the experimental data on which the hypothesis has been formulated.
furthermore, the careful analysis of counter examples found using this approach can lead to new hypotheses on how to reconcile the model with the experimental data.
we illustrate the usefulness of this approach by an example showing that we can infer a known interaction from the analysis of a model in which this interaction was purposefully ignored.
it is known that in mammalian keratinocytes, notch activates the expression of its ligand jagged.
for the purpose of this example, we ignore that information, and consider that notch inhibits the expression of the delta ligand, as it is the case in drosophila when building our original model.
such a model is not able to reproduce the wild type specification , and we have thus identified a gap between the hypothesized model and the experimental data.
the counter example traces we obtain indicate that the model is unable to reproduce the sustained notch signaling that is necessary in order to reproduce the observed wild type behavior.
this strongly suggests that the type of interaction between notch and its ligand in keratinocytes needs to be changed from inhibition to activation.
applying this change to the model and then running the verification algorithm again shows that the modified model is consistent with all specifications.
this means that our refinement process allows us to infer the hypothesis that was left out in this example.
furthermore, a similar approach, although using a different formalism, has lead to new biological insights about the vulval development in c. elegans  <cit> .
we believe that the possibility of in-silico execution of putative models is an efficient way to help biologists understand the behavior of a complex signaling networks.
counter example traces are helpful not only to detect that a model does not behave as expected, but also to provide a concrete basis on which one can reason about why this is indeed the case.
in addition, by querying the model, interesting executions could suggest new experiments to be tested in the lab.
we therefore expect that the application of our modeling methodology to other biological systems will help provide new biological insights.
an area in which our approach appears to be particularly useful is the analysis of cell-cell communication in multicellular models.
given the size of such models, efficient algorithms are required.
the efficient use of the modularity in multicellular models by our symbolic algorithm makes it particularly suitable for analyzing cell-cell communication processes together with complex intracellular pathways, a level of detail that would be intractable for enumerative or simpler symbolic algorithms.
the utility of this algorithm extends beyond the scope of our iterative improvement approach.
the computation of attractors of boolean networks is, for example, mainly done using enumerative methods.
such methods become intractable even for small models, and since boolean networks are a particular instance of qualitative networks, our symbolic algorithm can be directly applied to this problem.
modeling approach
in addition to being well suited for our formal verification based iterative improvement approach, using qualitative network has several benefits compared to using other existing frameworks, such a quantitative models using differential equations and boolean networks.
in this section, we discuss these advantages, and then stress the benefit of clearly separating the modeling of a system using a well defined framework like qualitative networks from its implementation.
using discrete qualitative variables sets, allows modeling at a level of abstraction which is similar to the one used in many biological studies.
quantitative methods, such as differential equation-based models, allow continuous and precise modeling of the execution of a system, but require fine-tuning of multiple parameters that are often not easily available for many biological systems.
the amount of precision provided by such methods may be unnecessary to analyze a putative model on the basis of experimental results that only provide qualitative levels for the components of the system.
the results obtained using qualitative network models of network motifs that occur frequently, such as the i1-ffl, indicate that the qualitative networks approach provides a realistic qualitative approximation of the experimentally observed behavior of the system.
the iterative improvement approach used for the construction of a consistent model as well as our analysis algorithms can be applied to boolean networks.
however, extending boolean networks allows both to create more detailed models in terms of granularity, and to model a richer set of biological interactions thanks to the introduction of the target function.
using more than two levels is necessary if an experimental result cannot be formulated as a requirement over boolean values.
this is the case if, for example, the knockout of gene a leads to an increase in the level of protein b, but the combined knockout of genes a and b leads to a significantly larger increase.
an additional example for which at least three levels are needed is the case when in one experiment the level of a protein decreases, while in another it increases.
the advantages of using qualitative networks are not limited to the increased expressive power of the specifications.
boolean networks are not sufficient to model some commonly observed patterns.
for example, the action of two transcription activators of a gene can be cumulative, whereas boolean networks only allow modeling the case in which the presence of one transcription factor is sufficient for transcription to happen.
furthermore, while the influence of the i1-ffl motif on the response time can be modeled using qualitative networks, this is not the case with boolean networks.
in the case of the studied example, the requirement that the presence of crp necessarily leads to the activation of galetk implies that, in a boolean network, the value of galetk will necessarily depend directly on the value of crp, and the contribution of gals needs to be ignored.
an even simpler example that cannot be reproduced by boolean networks is negative auto-regulation, where a transcription factor represses its own expression.
negative auto-regulation has been shown to lead to faster response times  <cit> .
in boolean networks, it is necessary to ignore negative auto regulation, since considering an inhibitory interaction of a component with itself would lead to oscillations.
these examples illustrate the need for extending boolean networks to the qualitative domain.
depending on the structure of the studied network and the corresponding specifications, boolean networks may however be sufficient in some cases.
since it is still possible to use the exact same analysis approach, boolean values should be used in such a situation to reduce the size of the state space.
possible extensions of the qualitative network framework include the addition of probabilities, resulting in markov processes, for which formal verification methods exist.
however, the necessary data for such models is often not available, and verification of probabilistic systems is computationally harder.
separation between the model and the implementation
we specify the molecular interactions between components by defining a qualitative network of the studied system.
compared to directly implementing an executable model, this approach allows to clearly separate this specification from the implementation.
this allows discussion at the level of the qualitative network, and thus leads to a improvement process at the level of abstraction often observed in diagrammatic models, which are commonly used by biologists.
furthermore, the clarity and coherence of the model is increased by defining a small set of target function that are consistently reused to model the same kind of interactions.
the benefits of a clear separation between the model and its implementation extend beyond qualitative networks.
using a modular approach based on simple building blocks allows a faster development of coherent executable models.
each type of building block represents a specific kind of biological interaction between components.
this approach separates the executable model into two distinct parts: the definition of the interactions between proteins such as activation or inhibition, and the internal implementation of the building blocks.
this modular design is particularly useful for fast search of multiple variations in the interaction model.
once the set of basic building blocks is implemented, interactions between biological elements of the model can be changed by replacing one building block by another, and new elements can be added as instances of the appropriate kind of building blocks, all this with minimal changes in the implementation of the model.
changes to the implementation of building blocks only need to be made in one place and can then be reused by all variations of the model.
in the case that we want to explore the behavior of the model when we use different ways for representing the interactions between biological elements, it is sufficient to replace one set of basic building blocks by another set of basic building blocks.
the clear separation between the interaction model and the implementation is also helpful for assessing the validity of the qualitative networks approach to modeling biological systems.
rather than having to understand the behavior of every single protein of every variation of the model, our approach allows to first agree on the definition of specific target functions, and then assess the plausibility of a particular qualitative network, whose representation is close to the common visualization of biochemical pathways.
the clear distinction between what we model and how we do it therefore makes the whole approach easier to understand.
the clear and concise definition of the behavior of a small set of basic building blocks makes it possible to evaluate  the conceptual elements of our approach.
finally, this approach allows to clearly separate methodological issues from issues that are only related to a particular qualitative network.
notions of concurrency
qualitative networks, like boolean networks, are synchronous.
at each time step, the state of all variables is updated according to the previous state of the system.
in this section, we discuss the applicability of other notions of concurrency in the context of qualitative networks.
the variables we use are qualitative representations of the level of protein activity or gene expression.
hence, these variables represent a large population of individual molecules.
the interactions between individual molecules are highly stochastic.
when using discrete levels of protein activity or gene expression, we consider that the biochemical reactions have a delay.
the stochasticity of the individual reactions results in small variations of this delay.
the level of all variables is updated according to the previous values.
a change therefore takes one time step to propagate from one variable to the next.
the time step can thus be considered as being the delay of the biochemical reactions.
in case there are major differences between the delays of the different reactions, it is possible to use target functions on the previous states of the systems rather than on the current state.
the impact of the small variations of the delay depend on the granularity of the model.
the minimal difference between two possible values of a variable is inversely proportional to the granularity of the variable.
in a model with a small granularity, variations of the delay would lead to the same value due to rounding.
in this case, the synchronous execution is sufficient to reproduce the behavior of the model.
as the granularity increases, rounding is no longer sufficient for handling the variability of the delays.
in order to model this variability, it is necessary to introduce non-determinism by allowing asynchronous execution of the model.
in a totally asynchronous system, a scheduler would choose any subset of the system which would be updated, while the other remain constant.
this means that the delay of an interaction is completely independent of the delay of other interactions of the system as well as of the previous delays of the same interaction.
this independence is thus an over-approximation of the expected behavior of the system, and is therefore likely to result in unrealistic executions.
therefore, while a synchronous execution is appropriate when the granularity is low, more precise models in terms of granularity call for a better handling of concurrency.
the state explosion problem
the exhaustive exploration of all possible executions of a system is exponential in the number of variables of the system.
this issue is called the state explosion problem.
for example, in the case of finding all infinitely visited states, the unrestricted number of initial states of a system with k variables and granularity n is nk.
enumerative methods are practically intractable even for small values of k. furthermore, the complexity of the notch/wnt model also caused the direct application of symbolic methods to fail for models of more than three cells.
the number of variables of a system not only depends on the size of the studied pathways, but also on the number of cells represented in the model.
hence, the number of initial states is also exponential in the number of cells, which makes multicellular models hard to analyze.
in order to be able to analyze larger models, it is therefore necessary to use more elaborate methods that have the potential for circumventing the state explosion problem in specific cases.
we provide an illustration of this approach in the symbolic algorithm used for computing the infinitely visited states.
we use both the modularity and the specific characteristics of transitions of our model in order to allow interleaving composition and computation of infinitely visited states.
the characteristics of our model also make it well suited for the use of partitioned transition relations.
we are able to compute the set of infinitely visited states for a  <dig> cells model, which has  <dig> ≈  <dig> initial states.
our approach is therefore well suited for analyzing multicellular qualitative networks.
this kind of model is particularly useful for studying pathways that involve inter-cellular communication.
in addition to using the specificity of the model for tailoring the verification algorithm, methods used in software and hardware modification, such as assume-guarantee reasoning  <cit> , or counterexample-guided abstraction refinement  <cit>  could also be applied to models of biological systems.
doing so has the potential for making it possible to analyze larger models.
in the context of the notch/wnt model, this could allow studying multidimensional tissues with larger numbers of cells, and thus the analysis of local effects, such as cell clusters.
CONCLUSIONS
in this work we propose an extension to boolean networks by allowing variables that represent the level of protein activity or gene expression to range over larger domains, combined with more flexible update options for these values.
we show that a qualitative networks model of a frequent network motif offers a valid qualitative approximation of the experimentally observed behavior of this motif.
this framework can be used to model recurrent networks.
in addition we use formal verification methods to analyze the steady state behavior of qualitative networks.
reasoning about sets of states rather than individual states, and using a partition reduction allows us to scale-up to very large models.
in particular, qualitative networks can be used to analyze multi-cellular models, and thus to study pathways that involve cell-cell communication.
similar to boolean networks, this approach could be useful for the analysis of biological pathways even where quantitative data is missing.
we believe that the ability to formally verify a biological model versus the laboratory observations offers great advantages in improving working hypotheses and suggesting new experimental directions that are likely to yield new findings.

 the emerging field of systems biology aims to gain high-level understanding of complex biological systems by studying the structure and dynamics of cellular functions. the construction and analysis of models describing biological processes are at the core of this new approach. computer science provides methods for simulating and analyzing these models. this allows performing in silico experiments that lead to new predictions and thus complement traditional experimental approaches. considering biology from a systems point of view also allows the use of methods designed for the engineering of artificial systems. a wide array of formalisms originating from computer science have already been used to model biological systems. these include process calculi, statecharts and live sequence charts, variants of petri nets, and hybrid models. recently, the term executable biology has been proposed for models that consist of an algorithm which mimics the behavior of a studied biological system. executable models of biological systems can be analyzed using formal verification methods, which were originally designed for the validation of computer systems. two possible approaches are testing and model checking. model checking consists of an exhaustive exploration of the state-space of a given system in order to verify that it always adheres to a set of requirements. in the case of biological systems, requirements are derived from experimental observations. the consistency of a putative executable model can thus be tested with respect to the data on which the model is based. in order to construct executable models of biological pathways, it is necessary to precisely define the interactions between biochemical components. one of the widely used formalism for representing interactions in biology are boolean networks. a boolean network consists of a set of boolean variables representing the activity of biochemical components such as proteins or genes. a boolean function is associated to every component. this function updates the value of the corresponding component at every step of the execution of a boolean network. in the context of biological pathways, a boolean network can be built from a graph representing biochemical interactions. every node of this graph represents a biochemical component. the arcs of the graph are used to represent the interactions between components. their weight can be positive or negative. this graph is used to define boolean functions for every component. the function sets the value of a component to true if the amount of activation is higher than the amount of inhibition. else the value of the component is set to false. boolean networks therefore provide a highly abstracted view of protein activity or gene expression. this level of abstraction often corresponds to the available qualitative experimental data. furthermore, the interactions between components are represented in a very similar way in the diagrammatic models that are widely used by biologists. boolean networks therefore appear to be well suited for the qualitative analysis of biological processes. they have been used recently for modeling a wide array of biological systems, such as the arabidopsis thaliana flower morphogenesis, the influence of physical stimuli on cell fate, the segment polarity gene network in drosophila melanogaster and the yeast cell cycle. boolean networks have been extended in several ways, including probabilistic boolean networks and the so-called generalized logical approach using multivalued levels of concentration. the generalized logical approach has been combined with formal verification methods, and constraint programming was used for finding steady states of large networks. in this work, we propose an extension of boolean networks called qualitative networks. we use this framework to construct an executable model which is analyzed using formal verification methods. we introduce two changes to the boolean networks framework. first, we use discrete variables instead of boolean variables for representing gene expression and protein activity levels. second, we capture a richer set of interactions between components than only activation and inhibition. any function of the current state of the system can be used to define the target towards which a component will move at the next step. at each step, the level of a component changes towards this target by at most one level. we use simulation and model checking to verify that the implementation of the qualitative network adheres to a set of specifications derived from the laboratory experimental observations. if necessary, we can incorporate new hypotheses into the model. once these are tested in-silico, they should also be validated experimentally. we propose that the stable states of the biological system correspond to the infinitely visited states of the executable model. we introduce an efficient symbolic algorithm for verifying if all infinitely visited states of a qualitative network satisfy the specifications. we consider the performance of this algorithm for models of different size and complexity, and show that it can be used to analyze large and complex qualitative networks. we illustrate the ability of qualitative networks to provide a valid qualitative representation of a studied biological system using a model of the incoherent type- feed-forward loop in the e. coli gal system. our approach to the modeling of signaling pathways allows the detection of gaps in the mechanistic understanding of the studied process, and suggests new biochemical interactions at a similar level of abstraction to the one observed in diagrammatic models. we illustrate this approach by a multicellular model of the interaction between the notch and the wnt pathways in mammalian skin cells.  qualitative networks in this study we propose an approach for constructing executable models of signaling pathways which is based on qualitative networks, an extension of boolean networks. in this section, we first explain how biological pathways can be modeled using boolean networks and then define qualitative networks. interactions in biological pathways can be represented as an edge-weighted interaction graph g. a node vi ∈ v represents a biological component such as a gene or a protein. an edge eij ∈ e from node vi to node vj has a weight αij, corresponding to the effect of the component represented by vi on the component represented by vj. activation is specified by a positive value of αij, inhibition by a negative value. this graph therefore represents the type of interaction between components as well as their strength in a similar way to the diagrammatic models used by biologists. li et al. propose the use of boolean networks for studying the interactions defined in such a graph. a boolean network b consists of a set of components c and a list of boolean functions f. each component ci ∈ c has a state, which can take two values: ci = if the component is active and ci = if the component is absent or inactive. we use k to represent the number of components. the state s of a boolean network corresponds to the vector ∈ { 1}k of the states of all components of the boolean network. a boolean function fi ∈ f is assigned to every component. this function maps the current state of the system to a boolean value representing the next state of the corresponding component. time in boolean networks is represented by a discrete variable t. we use the notations ci and s to represent the state of a component, respectively the state of the whole system, at time t. the value of each protein at the next time step is computed by applying the corresponding boolean function to the current state of the system: ci = fi )     given an interaction graph g, a boolean network b can be obtained as follows: each node of v corresponds to a component in c. the boolean function should reproduce the following behavior for every component: if the cumulative activation is higher than the cumulative inhibition, the state of the component is set to active and, vice-versa, if the inhibition is higher than the activation, it is set to inactive. hence, the boolean function corresponding to each component is defined as follows: fi={0if∑cj∈sαjicj<01if∑cj∈sαjicj>0ciif∑cj∈sαjicj=0 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgmbgzdawgaawcbagaemyaakgabeaakiabcicaoiabdohazjabcmcapiabg2da9maaceqabaqbaeqabmgaaaqaaiabicdawaqaaiabbmgapjabbagamnaaqafabaaccigae8xsde2aasbaasqaaiabdqgaqjabdmgapbqabagccqwgjbwydawgaawcbagaemoaaogabeaaaeaacqwgjbwydawgaaadbagaemoaaogabeaaliabgigiolabdohazbqab0gaeyyeiuoakiabgyda8iabicdawaqaaiabigdaxaqaaiabbmgapjabbagamnaaqafabagae8xsde2aasbaasqaaiabdqgaqjabdmgapbqabagccqwgjbwydawgaawcbagaemoaaogabeaaaeaacqwgjbwydawgaaadbagaemoaaogabeaaliabgigiolabdohazbqab0gaeyyeiuoakiabg6da+iabicdawaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaagcbagaeeyaakmaeeozay2aaabuaeaacqwfxoqydawgaawcbagaemoaaomaemyaakgabeaakiabdogajnaabaaaleaacqwgqbgaaeqaaaqaaiabdogajnaabaaameaacqwgqbgaaeqaasgaeyici4saem4camhabeqdcqghris5aogaeyypa0jaegimaadaaagaay5eaaaaaa@759c@ we extend the boolean network framework by using discrete domains to represent the state of a component and by allowing a rich set of interaction between components. the state of an active component can take several values, rather than being simply on. these values correspond to qualitative levels of gene expression or protein activity, for example, low, medium, high. this representation therefore is a middle ground between the boolean abstraction and quantitative models. the use of discrete domains requires to modify the way the updates of the network are defined. we introduce a target function that computes the discrete level towards which a component moves in the next step. the target function is not restricted to modeling activation and inhibition. every function of the current state of the system can be used as target function. this makes it possible to model interactions that cannot be represented in interaction graphs. one such example is complexation, where the amount of produced protein is bounded by the availability of the proteins required for the interaction. this situation can be modeled by using the minimum of the level of these proteins as a target function. more elaborate target functions can also be used to represent situations in which the effect of a specific activating or inhibiting edge is dependent on the level of additional proteins. we call this extended framework qualitative network. in the following paragraphs, we give a detailed description of qualitative networks. a qualitative network q consists of a set of components c and a list of target functions t. the state of a component ci ∈ c is an integer value between and n. these values represent qualitative levels of the components with the lowest and n the highest possible level for that component. for each component ci, there exists a target function targeti ∈ t which specifies the level towards which the value of the component should move. this function is computed on the state of the whole system. at each step, the value of the state of a component can only increase or decrease by a single level. the value of each component at the next time step is therefore computed as follows: ci={ci−1if targeti)<cici+1if targeti)>ciciif targeti)=ci      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabgucariabigdaxiabcmcapiabg2da9maaceqabaqbaeaabmgaaaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaogaeiikagiaemidaqnaeiykakiaeyoei0iaegymaedabagaeeyaakmaeeozaymaeeiiaaiaemidaqhcbigae8xyaemae8ncainae83zacmaemyzaumaemidaq3aasbaasqaaiabdmgapbqabagccqggoaakcqwgtbwucqggoaakcqwg0badcqggpaqkcqggpaqkcqgh8aapcqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabcmcapaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaogaeiikagiaemidaqnaeiykakiaey4kasiaegymaedabagaeeyaakmaeeozaymaeeiiaaiaemidaqnae8xyaemae8ncainae83zacmaemyzaumaemidaq3aasbaasqaaiabdmgapbqabagccqggoaakcqwgtbwucqggoaakcqwg0badcqggpaqkcqggpaqkcqgh+agpcqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabcmcapaqaaiabdogajnaabaaaleaacqwgpbqaaeqaaogaeiikagiaemidaqnaeiykakcabagaeeyaakmaeeozaymaeeiiaaiaemidaqnae8xyaemae8ncainae83zacmaemyzaumaemidaq3aasbaasqaaiabdmgapbqabagccqggoaakcqwgtbwucqggoaakcqwg0badcqggpaqkcqggpaqkcqgh9aqpcqwgjbwydawgaawcbagaemyaakgabeaakiabcicaoiabdsha0jabcmcapaaaaiaawuhaaiaaxmaacawljawaaewaaeaacqaiyagmaiaawicacaglpaaaaaa@9f11@ when n = this update function is equivalent to the update function of boolean networks. it follows that qualitative networks are a generalization of boolean networks: the qualitative network q is equivalent to the boolean network b. a particular kind of target function is directly obtained from the qualitative graph describing the studied biological system. it is used to model activation and inhibition in a similar way to the one describe for boolean networks. in boolean networks, the sum of all contributions was sufficient to decide on the next value of a component. here we need to obtain a target function that ranges over {..., n}, and for which intermediate values are meaningful. in order to obtain the target function for component ci, we compute separately the amount of activation acti and the amount of inhibition inhi sensed by this component. both of these values are scaled so that they are proportional to the weighted amount of activation, respectively inhibition. when all activating components are at their highest value, then acti = n. similarly, when all inhibiting components are at their highest value then inhi = n. acti=∑αji>0αjicj∑αji>0αjiinhi=∑αji<0αjicj∑αji<0αji mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakqaabeqaaiabdggahjabdogajjabdsha0naabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakiaeyypa0zaasaaaeaadaaeqbqaaggaciab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaogaem4yam2aasbaasqaaiabdqgaqbqabaaabagae8xsde2aasbaawqaaiabdqgaqjabdmgapbqabawccqgh+agpcqaiwaamaeqaniabgghildaakeaadaaeqbqaaiab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaaqaaiab=f7ahnaabaaameaacqwgqbgacqwgpbqaaeqaasgaeyopa4jaegimaadabeqdcqghris5aaaaaoqaaiabdmgapjabd6gaujabdigaonaabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakiaeyypa0zaasaaaeaadaaeqbqaaiab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaogaem4yam2aasbaasqaaiabdqgaqbqabaaabagae8xsde2aasbaawqaaiabdqgaqjabdmgapbqabawccqgh8aapcqaiwaamaeqaniabgghildaakeaadaaeqbqaaiab=f7ahnaabaaaleaacqwgqbgacqwgpbqaaeqaaaqaaiab=f7ahnaabaaameaacqwgqbgacqwgpbqaaeqaasgaeyipawjaegimaadabeqdcqghris5aaaaaaaa@7aa5@ the target function is the difference between the amount of activation and the amount of inhibition. when there is no inhibition, the output of the target function corresponds to the amount of activation. the level of the component will not exceed the amount of the activation. when the amount of inhibition is higher than the amount of activation, the value of the target function is zero. we need to consider a special case when all incoming edges of a component are inhibitions. in this case, we must assume that the mechanism leading to the presence of this component is not modeled, and that, in the absence of inhibition, the target function of this variable will return the maximal value. we obtain the following definition of the target function: targeti={max⁡−inhi)if max>0n−inhiif max≤0      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwg0badiegacqwfhbqycqwfybgccqwfnbwzcqwglbqzcqwg0baddawgaawcbagaemyaakgabeaakiabcicaoiabdohazjabcmcapiabg2da9maaceqabaqbaeaabigaaaqaaigbc2gatjabcggahjabciha4jabcicaoiabicdawiabcycasiabdggahjabdogajjabdsha0naabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakiaeyoei0iaemyaakmaemoba4maemiaag2aasbaasqaaiabdmgapbqabagccqggoaakcqwgzbwccqggpaqkcqggpaqkaeaacqqgpbqacqqgmbgzcqqggaaicqwftbqbcqwfhbqycqwf4baecqggoaakiigacqgfxoqydawgaawcbagaemoaaomaemyaakgabeaakiabcmcapiabg6da+iabicdawaqaaiabd6eaojabgkhitiabdmgapjabd6gaujabdigaonaabaaaleaacqwgpbqaaeqaaogaeiikagiaem4camnaeiykakcabagaeeyaakmaeeozaymaeeiiaaiae8xba0mae8xyaemae8heagnaeiikagiae4xsde2aasbaasqaaiabdqgaqjabdmgapbqabagccqggpaqkcqghkjyocqaiwaamaaaacagl7baacawljagaaczcamaabmaabagaeg4mamdacagloagaayzkaaaaaa@84bd@ we use qualitative networks to define executable models of biological systems. the list of target functions specifies how the different molecules in the system influence each other. when creating a model based on the mechanistic understanding of a biological process, we use this information to define a set of target functions. where necessary, we separate the representation of protein activity from the expression level of the corresponding gene. this is useful for modeling cases in which there are interactions at the transcription level as well as protein-protein interactions. during the analysis process, we may modify the definition of the target functions. while qualitative networks allow for any kind of target function, we use a subset of functions that represent specific biological interactions, such as the activation/inhibition function described above. this ensures that new hypothesized interactions between components are defined at a level where they have a clear biological meaning. the implementation of models defined as qualitative networks is described in the methods section. iterative model improvement a model in general, and an executable model in particular, is meant to represent how the system that we study actually works. thus, it is natural to expect that the model would be consistent with the data obtained from laboratory experimentation with this system. we use simulation and model checking, which allows to formally check all the different executions of the model against a formal specification. by exploring all the possible states and transitions of the model we can determine whether some requirement holds for the model. in the case that the property does not hold, the model-checking algorithm supplies a counter example: an execution of the system that does not satisfy a given requirement. we use the information provided by this counter example to suggest modifications of the qualitative network representing the studied system. our methodology in using qualitative networks is the iterative improvement process. we stress that the usage of qualitative networks is not restricted to this methodology. the advantages of multiple values, flexible target functions, and the analysis algorithms suggested below, can be used regardless of the modeling methodology chosen by the user. the iterative improvement process consists of the following. we first derive a set of specifications from prior laboratory experimental results. we build an initial putative model based on the current mechanistic understanding of the studied biological process and the working hypotheses we would like to evaluate. we then use formal verification methods to verify if the proposed model conforms with the specification. if the model fails this test, we deduce that the model does not conform to the data. we use counter examples, to suggest a revised model by making new hypotheses. the new model is then evaluated again, and these steps are repeated until no contradiction is found between the model and the specifications. at this point, we might have had to modify our hypotheses, both in discarding unverifiable assumptions and in adding assumptions that are necessary for the model to comply with the specifications. the resulting model, and the corresponding new set of hypotheses, are then known to be consistent with the experimental results. at this stage, it is possible to query the model to gain additional knowledge about the biological system. this can be done either by modifying the model to evaluate or by adding specifications that could provide additional information on the behavior of the system. hypotheses and changes in the model leading to interesting behaviors are candidates for experimental validation. the result of such an experiment, whether confirming the information gained from the computational analysis or not, can then be added to the specification. this creates a second, outer, improvement loop which combines computational and experimental validation. the specifications derived from laboratory experimental results can be separated into two categories. low-level specifications involve changes in the level of protein activity or gene expression. high-level specifications consider changes at the cellular level, such as cell fate or susceptibility to carcinogenesis. since we can only observe the level of individual proteins, we need to define how these relate to changes at the cellular level. experimental results are often related to manipulation, such as gene over-expression or knockout. these manipulation lead to changes in the behavior of the model. the methods section provides additional details on how specifications are defined and verified. this iterative improvement approach can be applied to highly non-deterministic models, hence the use of model checking. qualitative networks, like boolean networks, have a deterministic behavior, and it is therefore sufficient to perform one execution per initial state, and verify if each of these executions satisfies the specifications. if some executions do satisfy the specifications, but others do not, then the model is not constrained enough, and we need to formulate hypotheses that allow avoiding the executions which do not satisfy the specifications. ideally, the new hypotheses would not contradict the current mechanistic understanding of the system, but rather suggest additional putative interactions between components. if there is a specification which is not satisfied by any execution of the model, then there is a contradiction between the mechanistic understanding of the system and the experimental data. in this case, the modifications needed for the model to satisfy the requirement will suggest alternative mechanistic models of the studied system that are more consistent with the experimental data. in the following sections, we describe two additional methods for analyzing the model during the improvement process: the use of non-determinism on parts of the model, and the efficient computation of all infinitely visited states. use of non-determinism in the case where every interaction in the model is precisely defined, enumerative networks are deterministic. non-determinism can, however, be used to represent unknown interactions. in particular, rather than studying a complete multicellular model, the model of a single cell can be studied by representing inter-cellular interactions non-deterministically. the level of proteins that are external to the cell are modeled by non-deterministic variables. these variables cannot jump arbitrarily from one value to another value, but only change by one level per time step. at each time step, the variable can either increase, decrease or remain constant. this is sufficient to capture any possible behavior since updates of components are bounded by equation the set of possible behaviors of the state of this protein is therefore a superset of the possible behaviors for any possible target function. instead of simulation, we use model checking, which explores all possible executions of a non-deterministic system. this leads to the analysis of a superset of the possible executions of any model with more constrained inter-cellular interactions. if no execution of the system satisfies the specification, then it is impossible that any model including this cell can satisfy the specifications, no matter how the inter-cellular interactions are defined. if some executions satisfy the specifications, but other do not, then it is not possible to predict the behavior of the cell with more constrained inter-cellular interactions. this approach is therefore useful for finding out if the improvement should be performed at the level of the inter-cellular communication or at the intracellular level. finding infinitely visited states we propose an efficient symbolic algorithm for computing the infinitely visited states of a qualitative network. a state s is said to be infinitely visited if there is an execution of the network in which the state appears infinitely often. in a biological context, the infinitely visited states correspond to the stable states of the system. in the context of specification-based analysis of the model, we are therefore interested in knowing if the specification holds on all infinitely visited states. this captures the fact that, as opposed to computerized systems, biological system do not have an initial state, but rather are continuously behaving. after cell division, for example, the two resulting cells would have proteins levels that are at a qualitative level similar to the levels in the original cell. states that are not infinitely visited can thus be considered as unstable states. if the model is in such a state, it will evolve towards an infinitely visited state in a bounded number of steps, and then remain in infinitely visited states for an indefinite time. therefore, we do not insist that the specification holds on states that are not infinitely visited. analysis of boolean networks also include studies of the attractors of the model. in discrete, deterministic models, attractors are loops of one or more states that are visited infinitely often. questions of interest include: does the system end up in one or several attractors, and what initial conditions lead to which attractor. this analysis requires exploring the executions starting from all initial states. the number of possible initial states is exponential in the number of variables, and therefore we try to prune them by using the available information on the system. in case that no biologically meaningful definition of the set of initial state exists, it is necessary to consider all states as being initial. in this case, enumerative exploration is practically intractable even for relatively small models. in the methods section, we describe a symbolic algorithm for finding all the attractors of a model, based on compact representation of sets of states using binary decision diagrams. the algorithm we propose uses the structure of multicellular qualitative networks by interleaving composition and computation of infinitely visited states. we first consider a partial model consisting only of the first cell. the behavior of components in neighbor cells are abstracted by allowing them to change non-deterministically. we compute the set of infinitely visited states of this partial model. this set contains the projection of any infinitely visited state of the complete model onto the variables of the partial model. we use this information to compute the infinitely visited states of a partial model consisting of the first two cells. this process is repeated to obtain the complete model and the set of infinitely visited states. we further optimize this algorithm by using partitioned transition relations. based on the target functions list of the qualitative network, we define a transition relation for a single cell. we then use it repeatedly to compute the next state of the multicellular system, rather than using a single transition relation for the whole model. the methods section contains a detailed description of this algorithm. modeling network motifs using qualitative networks transcription regulation networks contain patterns that occur significantly more often in a real network than in a random network with the same characteristics. these patterns are called network motifs. the multiple values used by qualitative networks may be required to accurately model the qualitative behavior of common network motifs in a way that would be difficult, if not impossible, to achieve with boolean values. we illustrate this ability by considering the case of the incoherent type- feed-forward loop in e. coli. in this motif, transcription factor x activates transcription factor y and gene z, and transcription factor y represses z. the response time is defined as the time needed for z to increase halfway to its steady-state level. compared to simple regulation, the i1-ffl has a shorter response time to the activation of x. when x is activated, the level of both y and z start to increase. once the level of y exceeds the repression threshold for the z promoter, it will start repressing z. as a consequence, the level of z starts decreasing towards a low steady-state level. the response time is thus shorter, because the steady-state level of z is lower than for simple activation, but its initial transcription rate is the same in both cases. this behavior has been observed both using quantitative models and experimentally in the gal system in e. coli. in the gal system, an i1-ffl links the activator crp to the galactose utilization operon galetk. when there is glucose starvation, crp gets activated. crp then activates both galetk and gals, and gals inhibits galetk, hence the i1-ffl motif. in addition, β-d-galactose or d-fucose unbind gals from the galetk promoter region, thus cancelling its inhibitory effect. finally, gals negatively auto-regulates its own expression. we construct a qualitative network model of this system based on these interactions. this model needs to satisfy two requirements to be realistic. first, the activation of crp must lead to the activation of galetk independently of the level of galsinh. second, the time needed by galetk to increase halfway to its steady-state needs to be shorter when galsinh is absent than when it is present. figure 2b shows the evolution over time of the level of galetk after crp is activated for both case. this result is qualitatively similar to the experimental measurements of the galetk expression performed by mangan et al.. qualitative networks are able to provide realistic qualitative approximations of the behavior of further common network motifs, such as other types of feed-forward loops and regulated-feedback motifs that appear in developmental transcription networks. case study: a model of the interaction between the notch and wnt signaling pathways we apply our modeling approach to a model representing the crosstalk between the notch and wnt pathways in mammalian keratinocytes. both pathways play a key role in the control of cell proliferation and differentiation and have been linked to several types of cancer. we first present the biological background, then introduce the qualitative network model for a single cell and for a model with five cell representing the multiple layers of the mammalian epidermis. we define a set of specifications derived from laboratory experimental perturbations on these pathways, and use them to analyze the possible executions of the model. crosstalk between the notch and wnt signaling pathways notch signaling plays a key role in the control of the development of various tissues, with the particularity that, depending on the context, it can either induce differentiation or maintain cells in a proliferating state. the notch receptor protein is located across the cell membrane and is composed of two parts – an extracellular component and a cytoplasmic part. the extracellular component of notch can bind to a ligand expressed on the membrane of a neighboring cell. in mammals, there are four types of notch proteins, three types of delta-like ligands and two types of serrate-like ligands, which are called jagged. the binding of a ligand to the notch receptor leads to a proteolytic cleavage of the two parts of the notch protein, thus releasing its cytoplasmic part – the notch intracellular domain. notch-ic then travels to the nucleus, where it binds to the transcription factor csl. in the absence of notch-ic, csl acts as a transcription inhibitor, but the presence of notch-ic converts it to a transcription activator, thus leading to the expression of its transcriptional targets. the wnt signaling pathway has a critical regulating role in stem cells, and has been associated with cancer in several tissues. the central player in the canonical wnt signaling pathway is the protein β-catenin. in the absence of wnt signaling, the newly formed β-catenin is found in a degradation complex composed of the proteins adenomous polyposis coli, the scaffolding protein axin and the glycogen-synthase kinase-3β. the gsk3β phosphorilates β-catenin, which is then degradated by proteasomes. wnt is a short range signaling protein which interacts with a serpentine receptor of the frizzled family. the consequence of this interaction are not completely understood, but it appears that an axin-binding molecule, dishevelled, gets phosphorylated and starts inhibiting the degradation of β-catenin. this β-catenin can therefore accumulate in the cytoplasm and travel to the nucleus, where it interacts with the tcf/lef dna-binding proteins. in the absence of β-catenin, a co-repressor like groucho binds to tcf/lef, thus blocking the transcription of its downstream target genes. when interacting with β-catenin, tcf/lef becomes a transcription activator, thus leading to the activation of the downstream target genes of the wnt signaling pathway. the outermost layer of the human skin, the epidermis, can be divided into several layers: the basal layer, which is the closest to the underlying dermis, the suprabasal layer, and the cornified layer, which is mainly composed of dead cells. in normal tissue, most proliferating keratinocytes are located in the basal layer. they migrate to the suprabasal layer when they have initiated terminal differentiation and then stop proliferating. notch is expressed in the human epidermis, at a higher level in the suprabasal layers than in the basal layer. the jagged ligand is co-expressed with notch in the skin and has been proposed as a positive feedback mechanism for sustained notch signaling triggering differentiation. the interaction between the notch and the wnt pathways are involved in the transition from a proliferating to a differentiated state in keratinocytes. in vitro experiments have shown that notch signaling induces terminal differentiation in both murine and human skin. the tissue specific ablation of the notchl receptor in mouse epidermis leads to hyperplasia in the epidermis, the development of tumors and an increased susceptibility to chemical carcinogens – thus indicating that notch acts as a tumor suppressor in mammal skin. notch was suggested to interact with the wnt pathway through its direct transcriptional target p which is a negative transcription regulator of wnt. model construction we build a model of the notch and wnt pathways in the epidermal layer of mammalian skin. our model is composed of five identical keratinocytes, each of them representing one layer of the epidermis. protein activity and gene expression are represented by variables that can take four values: off, low, medium and high. all activation and inhibition interactions have an equal weight. figure represents the two pathways in a single cell. except for the intracellular part of notch, the target function of all other components is the difference between activation and inhibition, as defined by equation the target function of the notch-ic component is the minimum between the level of notch receptor on the cell membrane and the level of ligand on the membrane of the neighboring cells. the complete model consists of a single row of five cells, which are named cell to cell neighboring cells are connected both through the wnt and the notch pathways. cell the leftmost cell on the visualization of the model, represents the lowest layer of keratinocytes of the epidermis, which means that the dermis is to its immediate left. this cell should adopt a proliferating fate. the upper skin, to which the keratinocytes migrate before they die, is represented on the right side of the model. we consider that the two immediate neighbors contribute to the level of ligand sensed by a cell. the level of ligand sensed by a cell is activated by the level of ligand on the membrane of its immediate neighbor, and is thus maximal when the level of ligand in both neighbors is high. the level of jagged on both extremities of the model is constant at medium. similarly, we connect the level of wnt produced by a cell to the level of wnt sensed by the nearby cells. we consider the case in which a cell does not sense wnt emitted by itself. in this case we can use the same double activation scheme as for the ligands. on the left side of the model, the level of wnt is fixed to high since the dermis is known to emit wnt signaling proteins. in contrast, the level on the right side is fixed to low since wnt signaling does not occur in the upper layers of the skin. specifications in order to be able to define high level specifications, we need to map changes in the level of protein activity to a certain cell fate. we consider the balance between the downstream target genes of wnt signaling and the downstream target genes of notch signaling. the downstream targets of the wnt signaling are known to maintain the cells in a proliferating state, whereas the downstream targets of notch initiate differentiation. therefore, we make the assumption that a cell in which gt > gt is proliferating, whereas a cell in which gt < gt is differentiated. furthermore, since an increase in the level of gt or a decrease in the level of gt makes a cell more likely to be proliferating, we assume that this change also makes the cell more sensitive to chemically-induced carcinogenesis. for experiments involving the knockout or the over-expression of a gene, we compare the level of components when the wild type model is in a stable state to the level of the proteins when the altered model reaches a stable state again. the first high level specification defines the wild type state, and can be split into two parts: cell is proliferating, and cell4– are differentiated. h: gt > gt in cell1 h: gt < gt in cell4–5 we do not use cell2– in these specification, but if a cell has adopted a differentiated fate, cells in a higher sub-layer cannot be proliferating. we formulate this as an additional requirement: h gt = gt in celli ⇒ ∀j > i, gt ≤ gt in cellj. h gt < gt in celli ⇒ ∀j > i, gt < gt in cellj. we derive further specifications from the experiments performed by nicolas et al. on mouse keratinocytes. they note that notch knockout leads to increased proliferation, which we translate into a requirement that cell is also proliferating. h3: notch = ko ⇒ gt > gt in cell4 after several months, mice in which notch is knocked out begin to develop basal cell-like carcinoma. they are also significantly more sensitive to chemically induced carcinogenesis. in order to relate these long-term changes to variations at the protein level, we assume that the short term effect of the notch knockout is an increase of gt in several cells, or a decrease of gt in several cells, or both. h4: notch = ko ⇒ gt increases in cell1– or gt decreases in cell1–5 the impact of p knockout was also considered. mice with p deficiency were more sensitive to chemically induced carcinogenesis, but did not develop tumors on their own. we relate this result to the level of gt and gt in a similar way to notch knockout. h5: p = ko ⇒ gt increases in cell1– or gt decreases in cell1–5 the work by nicolas et.al. on mice keratinocytes, as well as the work of devgan et.al. on the link between notch and wnt through p also leads to several low level specifications. these specifications were directly considered when constructing the model and are not detailed here. analysis of the model we compute the set of infinitely visited states of the model described above. we obtain a total of states. all of these states adhere to the specifications h and h these states, however, differ as far as the fate of cell is concerned. this cell is either proliferating, differentiated or in transition from proliferating to differentiated. several of these attractors consist of more than one state, and the relation between gt and gt in cell can change between the different states of a single attractor. we can therefore conclude that, in this model, cell is no longer receiving constant signals for proliferation, but has not yet committed to terminal differentiation. after the model reaches a stable state, modifying it to knockout notch or p allows to verify that the model also satisfies specifications h respectively h further insights can be gained by studying variations of the model. we consider the hypothesis that notch-ic activates the transcription of β-catenin, thus partially compensating its inhibitory influence on β-catenin through the wnt pathway. this hypothesis can be rejected by analyzing a single cell model with non-deterministic levels for both external wnt and sensed ligand. no execution of this model can satisfy requirement h for all the states of the execution that occur infinitely often. this means that in any multicellular model with this interaction, no cell could be proliferating. we can therefore conclude that some other mechanism activating β-catenin must exist. notch-ic may play a role in the activation of β-catenin, but cannot be the only activating protein. we also consider the hypothesis that notch signaling uses the delta ligand instead of the jagged ligand. while there is a positive feedback mechanism from notch-ic to jagged, notch-ic inhibits the expression of the delta ligand on the cell surface. we find that such a model cannot reproduce the wild type behavior. by analyzing the counterexample traces, we notice that due to the inhibition of the ligand, sustained notch signaling cannot be achieved in the higher layers of the epidermis. this is consistent with the experimental analysis of the level of delta in the epidermis showing that delta is expressed in the basal layer only. in their analysis, lowell et al. describe stem-cell clusters with local concentration of delta. in order to model and analyze local effect such as cell clusters, it would be necessary to consider significantly larger multidimensional models. performance analysis in this section, we describe the performance of the symbolic algorithm used to compute all infinitely visited states of a qualitative network. we first consider variations of the notch/wnt model described above. we then use an arbitrary model to show that the method can be applied to models of more complex interactions between components. finally, we show that our algorithm performs better than simpler symbolic algorithms. the methods section gives additional details about the settings used for these analyses. the computation of all infinitely visited states of the five-cells notch/wnt model takes minutes. this model has ≈ initial states. this is a significantly larger number than previously studied by boolean networks. in order to evaluate the performance of our symbolic algorithm on larger models, we extend the notch/wnt model up to cells. we add cells both to the basal layer and to the suprabasal layer. as in the five-cells model, the level of notch receptor is set to off in the basal layer, and to high in the suprabasal layer. while the original model has only one cell with an intermediate level of notch receptor, we introduce a second cell in the larger model in order to obtain a gradient. these models have been developed for the purpose of performance analysis only, and do not provide additional biological insights. table shows the time needed to compute the set of all infinitely visited states for models with to cells. the results illustrate that we can compute the set of infinitely visited states of a model which has more than initial states. it can be observed that adding a cell in the basal layer has a significantly smaller impact on the performance than adding a cell in the suprabasal layer. this is due to the fact that the possible behaviors of the cells in the basal layer are very limited due to the lack of notch signaling. in general, the performance of symbolic algorithms is dependent on how compact the representation of the sets of states and the transition relation are. this complexity is directly dependent on the structure of the model. in the case of qualitative networks, choosing a different combination of target functions has an impact on the performance even though the number of components stays the same. we illustrate this by constructing an arbitrary qualitative network, which has no biological meaning except that we only use the activation/inhibition target function. this model contains the same number of components per cell as the notch/wnt model. the interactions between the components are more complex. in order to compare the complexity of the two models, we consider the size of the bdd used to represent the transition relation of a single cell. this size, measured as the number of nodes in the bdd, relates more or less to the required memory. the arbitrary qualitative network is represented using 279' bdd nodes, compared to 39' bdd nodes for the notch/wnt model. this increase in complexity is reflected in the performance of the infinitely visited states computation. when considering only three cells this computation takes seconds for the notch/wnt model, and minutes for the arbitrary qualitative network. our algorithm interleaves the composition of a multi-cellular model with the computation of infinitely visited states. this approach outperforms symbolic algorithms that directly compute the set of infinitely visited states of the complete model. the application of the direct algorithm to models of the notch/wnt pathway with four or more cells fails, because the symbolic representations of the sets of states becomes too large. on the three-cell model, our algorithm needs seconds and the direct algorithm seconds. the advantage of first computing the infinitely visited states on partial models is that the sets of potentially infinitely visited states is smaller than in the direct algorithm. this is illustrated in table which provides details about the application of our method to the ten-cells notch/wnt model. for each partial model, we indicate the initial set of states we start from, as well as the set of infinitely visited states of the partial model after the computation. in particular, it can be noted that the computation of the infinitely visited states of the complete model starts from a set of approximately states. an algorithm that directly computes all infinitely visited states would start from the set of all states of the system. the symbolic approach we propose is able to compute the infinitely visited states for large multicellular qualitative networks and to handle complex models. by tailoring the algorithm specifically for qualitative networks, we can analyze larger models than by directly computing the infinitely visited states. discussion in this work, we extend boolean networks to qualitative networks; networks that allow variables to take discrete values representing qualitative levels of protein activity or gene expression. by introducing a target function, we allow a rich set of interactions between components. we propose a method for analyzing qualitative networks, that uses formal verification of specifications derived from experimental observations for iteratively improving the model. in order to efficiently analyze the properties of infinitely visited states of large models, we introduce a symbolic algorithm which builds on the structure of our modeling method. we apply this method to build and analyze a multicellular model of the interactions between the wnt and notch pathways in mammalian keratinocytes. we show that the proposed symbolic algorithm can be applied to large and complex multicellular models. model analysis approach our modeling approach is aimed at using formal verification methods for the analysis of signaling networks in an iterative improvement process which is based on specifications that are directly derived from experimental results. this work therefore combines formal verification methods with the widely used formalism of boolean network, and its extension to qualitative domains. since boolean networks are a particular case of qualitative networks, the iterative improvement approach can be directly applied to boolean networks as well. our method allows to build a model at a level of abstraction that is similar to the one commonly used in cell biology, and to analyze it using specifications that are directly derived from experimental studies. this approach is therefore useful for biologists to identify gaps between the hypothesized mechanistic description of a pathway of interest, and the experimental data on which the hypothesis has been formulated. furthermore, the careful analysis of counter examples found using this approach can lead to new hypotheses on how to reconcile the model with the experimental data. we illustrate the usefulness of this approach by an example showing that we can infer a known interaction from the analysis of a model in which this interaction was purposefully ignored. it is known that in mammalian keratinocytes, notch activates the expression of its ligand jagged. for the purpose of this example, we ignore that information, and consider that notch inhibits the expression of the delta ligand, as it is the case in drosophila when building our original model. such a model is not able to reproduce the wild type specification, and we have thus identified a gap between the hypothesized model and the experimental data. the counter example traces we obtain indicate that the model is unable to reproduce the sustained notch signaling that is necessary in order to reproduce the observed wild type behavior. this strongly suggests that the type of interaction between notch and its ligand in keratinocytes needs to be changed from inhibition to activation. applying this change to the model and then running the verification algorithm again shows that the modified model is consistent with all specifications. this means that our refinement process allows us to infer the hypothesis that was left out in this example. furthermore, a similar approach, although using a different formalism, has lead to new biological insights about the vulval development in c. elegans. we believe that the possibility of in-silico execution of putative models is an efficient way to help biologists understand the behavior of a complex signaling networks. counter example traces are helpful not only to detect that a model does not behave as expected, but also to provide a concrete basis on which one can reason about why this is indeed the case. in addition, by querying the model, interesting executions could suggest new experiments to be tested in the lab. we therefore expect that the application of our modeling methodology to other biological systems will help provide new biological insights. an area in which our approach appears to be particularly useful is the analysis of cell-cell communication in multicellular models. given the size of such models, efficient algorithms are required. the efficient use of the modularity in multicellular models by our symbolic algorithm makes it particularly suitable for analyzing cell-cell communication processes together with complex intracellular pathways, a level of detail that would be intractable for enumerative or simpler symbolic algorithms. the utility of this algorithm extends beyond the scope of our iterative improvement approach. the computation of attractors of boolean networks is, for example, mainly done using enumerative methods. such methods become intractable even for small models, and since boolean networks are a particular instance of qualitative networks, our symbolic algorithm can be directly applied to this problem. modeling approach in addition to being well suited for our formal verification based iterative improvement approach, using qualitative network has several benefits compared to using other existing frameworks, such a quantitative models using differential equations and boolean networks. in this section, we discuss these advantages, and then stress the benefit of clearly separating the modeling of a system using a well defined framework like qualitative networks from its implementation. using discrete qualitative variables sets, allows modeling at a level of abstraction which is similar to the one used in many biological studies. quantitative methods, such as differential equation-based models, allow continuous and precise modeling of the execution of a system, but require fine-tuning of multiple parameters that are often not easily available for many biological systems. the amount of precision provided by such methods may be unnecessary to analyze a putative model on the basis of experimental results that only provide qualitative levels for the components of the system. the results obtained using qualitative network models of network motifs that occur frequently, such as the i1-ffl, indicate that the qualitative networks approach provides a realistic qualitative approximation of the experimentally observed behavior of the system. the iterative improvement approach used for the construction of a consistent model as well as our analysis algorithms can be applied to boolean networks. however, extending boolean networks allows both to create more detailed models in terms of granularity, and to model a richer set of biological interactions thanks to the introduction of the target function. using more than two levels is necessary if an experimental result cannot be formulated as a requirement over boolean values. this is the case if, for example, the knockout of gene a leads to an increase in the level of protein b, but the combined knockout of genes a and b leads to a significantly larger increase. an additional example for which at least three levels are needed is the case when in one experiment the level of a protein decreases, while in another it increases. the advantages of using qualitative networks are not limited to the increased expressive power of the specifications. boolean networks are not sufficient to model some commonly observed patterns. for example, the action of two transcription activators of a gene can be cumulative, whereas boolean networks only allow modeling the case in which the presence of one transcription factor is sufficient for transcription to happen. furthermore, while the influence of the i1-ffl motif on the response time can be modeled using qualitative networks, this is not the case with boolean networks. in the case of the studied example, the requirement that the presence of crp necessarily leads to the activation of galetk implies that, in a boolean network, the value of galetk will necessarily depend directly on the value of crp, and the contribution of gals needs to be ignored. an even simpler example that cannot be reproduced by boolean networks is negative auto-regulation, where a transcription factor represses its own expression. negative auto-regulation has been shown to lead to faster response times. in boolean networks, it is necessary to ignore negative auto regulation, since considering an inhibitory interaction of a component with itself would lead to oscillations. these examples illustrate the need for extending boolean networks to the qualitative domain. depending on the structure of the studied network and the corresponding specifications, boolean networks may however be sufficient in some cases. since it is still possible to use the exact same analysis approach, boolean values should be used in such a situation to reduce the size of the state space. possible extensions of the qualitative network framework include the addition of probabilities, resulting in markov processes, for which formal verification methods exist. however, the necessary data for such models is often not available, and verification of probabilistic systems is computationally harder. separation between the model and the implementation we specify the molecular interactions between components by defining a qualitative network of the studied system. compared to directly implementing an executable model, this approach allows to clearly separate this specification from the implementation. this allows discussion at the level of the qualitative network, and thus leads to a improvement process at the level of abstraction often observed in diagrammatic models, which are commonly used by biologists. furthermore, the clarity and coherence of the model is increased by defining a small set of target function that are consistently reused to model the same kind of interactions. the benefits of a clear separation between the model and its implementation extend beyond qualitative networks. using a modular approach based on simple building blocks allows a faster development of coherent executable models. each type of building block represents a specific kind of biological interaction between components. this approach separates the executable model into two distinct parts: the definition of the interactions between proteins such as activation or inhibition, and the internal implementation of the building blocks. this modular design is particularly useful for fast search of multiple variations in the interaction model. once the set of basic building blocks is implemented, interactions between biological elements of the model can be changed by replacing one building block by another, and new elements can be added as instances of the appropriate kind of building blocks, all this with minimal changes in the implementation of the model. changes to the implementation of building blocks only need to be made in one place and can then be reused by all variations of the model. in the case that we want to explore the behavior of the model when we use different ways for representing the interactions between biological elements, it is sufficient to replace one set of basic building blocks by another set of basic building blocks. the clear separation between the interaction model and the implementation is also helpful for assessing the validity of the qualitative networks approach to modeling biological systems. rather than having to understand the behavior of every single protein of every variation of the model, our approach allows to first agree on the definition of specific target functions, and then assess the plausibility of a particular qualitative network, whose representation is close to the common visualization of biochemical pathways. the clear distinction between what we model and how we do it therefore makes the whole approach easier to understand. the clear and concise definition of the behavior of a small set of basic building blocks makes it possible to evaluate the conceptual elements of our approach. finally, this approach allows to clearly separate methodological issues from issues that are only related to a particular qualitative network. notions of concurrency qualitative networks, like boolean networks, are synchronous. at each time step, the state of all variables is updated according to the previous state of the system. in this section, we discuss the applicability of other notions of concurrency in the context of qualitative networks. the variables we use are qualitative representations of the level of protein activity or gene expression. hence, these variables represent a large population of individual molecules. the interactions between individual molecules are highly stochastic. when using discrete levels of protein activity or gene expression, we consider that the biochemical reactions have a delay. the stochasticity of the individual reactions results in small variations of this delay. the level of all variables is updated according to the previous values. a change therefore takes one time step to propagate from one variable to the next. the time step can thus be considered as being the delay of the biochemical reactions. in case there are major differences between the delays of the different reactions, it is possible to use target functions on the previous states of the systems rather than on the current state. the impact of the small variations of the delay depend on the granularity of the model. the minimal difference between two possible values of a variable is inversely proportional to the granularity of the variable. in a model with a small granularity, variations of the delay would lead to the same value due to rounding. in this case, the synchronous execution is sufficient to reproduce the behavior of the model. as the granularity increases, rounding is no longer sufficient for handling the variability of the delays. in order to model this variability, it is necessary to introduce non-determinism by allowing asynchronous execution of the model. in a totally asynchronous system, a scheduler would choose any subset of the system which would be updated, while the other remain constant. this means that the delay of an interaction is completely independent of the delay of other interactions of the system as well as of the previous delays of the same interaction. this independence is thus an over-approximation of the expected behavior of the system, and is therefore likely to result in unrealistic executions. therefore, while a synchronous execution is appropriate when the granularity is low, more precise models in terms of granularity call for a better handling of concurrency. the state explosion problem the exhaustive exploration of all possible executions of a system is exponential in the number of variables of the system. this issue is called the state explosion problem. for example, in the case of finding all infinitely visited states, the unrestricted number of initial states of a system with k variables and granularity n is nk. enumerative methods are practically intractable even for small values of k. furthermore, the complexity of the notch/wnt model also caused the direct application of symbolic methods to fail for models of more than three cells. the number of variables of a system not only depends on the size of the studied pathways, but also on the number of cells represented in the model. hence, the number of initial states is also exponential in the number of cells, which makes multicellular models hard to analyze. in order to be able to analyze larger models, it is therefore necessary to use more elaborate methods that have the potential for circumventing the state explosion problem in specific cases. we provide an illustration of this approach in the symbolic algorithm used for computing the infinitely visited states. we use both the modularity and the specific characteristics of transitions of our model in order to allow interleaving composition and computation of infinitely visited states. the characteristics of our model also make it well suited for the use of partitioned transition relations. we are able to compute the set of infinitely visited states for a cells model, which has ≈ initial states. our approach is therefore well suited for analyzing multicellular qualitative networks. this kind of model is particularly useful for studying pathways that involve inter-cellular communication. in addition to using the specificity of the model for tailoring the verification algorithm, methods used in software and hardware modification, such as assume-guarantee reasoning, or counterexample-guided abstraction refinement could also be applied to models of biological systems. doing so has the potential for making it possible to analyze larger models. in the context of the notch/wnt model, this could allow studying multidimensional tissues with larger numbers of cells, and thus the analysis of local effects, such as cell clusters.  in this work we propose an extension to boolean networks by allowing variables that represent the level of protein activity or gene expression to range over larger domains, combined with more flexible update options for these values. we show that a qualitative networks model of a frequent network motif offers a valid qualitative approximation of the experimentally observed behavior of this motif. this framework can be used to model recurrent networks. in addition we use formal verification methods to analyze the steady state behavior of qualitative networks. reasoning about sets of states rather than individual states, and using a partition reduction allows us to scale-up to very large models. in particular, qualitative networks can be used to analyze multi-cellular models, and thus to study pathways that involve cell-cell communication. similar to boolean networks, this approach could be useful for the analysis of biological pathways even where quantitative data is missing. we believe that the ability to formally verify a biological model versus the laboratory observations offers great advantages in improving working hypotheses and suggesting new experimental directions that are likely to yield new findings. 
we consider networks in which elements range over a small finite domain allowing more flexibility than boolean values, and add target functions that allow to model a rich set of behaviors. if the model does not conform to the data, we suggest a revised model and the new hypotheses are tested in-silico. here we introduce qualitative networks, an extension of boolean networks. we introduce an approach for constructing computational models of biological systems that extends the framework of boolean networks and uses formal verification methods for the analysis of the model. hence, this approach has the potential to efficiently complement experimental studies in biology. this approach can scale to multicellular models of complex pathways, and is therefore a useful tool for the analysis of complex biological systems. we propose a symbolic algorithm for analyzing the steady state of these networks, allowing us to scale up to a system consisting of elements and state spaces of approximately states. we illustrate the usefulness of this approach through a model of the interaction between the notch and the wnt signaling pathways in mammalian skin, and its extensive analysis. with this framework, we use formal verification methods to check whether a model is consistent with the laboratory experimental observations on which it is based.
 a central goal of systems biology is to model and analyze biological signaling pathways that interact with one another to form complex networks. here we introduce qualitative networks, an extension of boolean networks. with this framework, we use formal verification methods to check whether a model is consistent with the laboratory experimental observations on which it is based. if the model does not conform to the data, we suggest a revised model and the new hypotheses are tested in-silico.  we consider networks in which elements range over a small finite domain allowing more flexibility than boolean values, and add target functions that allow to model a rich set of behaviors. we propose a symbolic algorithm for analyzing the steady state of these networks, allowing us to scale up to a system consisting of elements and state spaces of approximately states. we illustrate the usefulness of this approach through a model of the interaction between the notch and the wnt signaling pathways in mammalian skin, and its extensive analysis.  we introduce an approach for constructing computational models of biological systems that extends the framework of boolean networks and uses formal verification methods for the analysis of the model. this approach can scale to multicellular models of complex pathways, and is therefore a useful tool for the analysis of complex biological systems. the hypotheses formulated during in-silico testing suggest new avenues to explore experimentally. hence, this approach has the potential to efficiently complement experimental studies in biology. 