Basins of Attraction, Commitment Sets and Phenotypes of Boolean Networks
Hannes Klarner Frederike Heinitz Sarah Nee Heike Siebert





Basins of Attraction, Commitment Sets and Phenotypes of Boolean Networks
    Hannes Klarner, Frederike Heinitz, Sarah Nee, Heike Siebert

Department of Mathematics, Freie Universität Berlin, Germany

E-mail: Hannes.Klarner@FU-Berlin.de
    December 30, 2023
==================================================================================================================================================================





The attractors of Boolean networks and their basins have been shown to be highly relevant for model validation and predictive modelling, e.g., in systems biology. 
Yet there are currently very few tools available that are able to compute and visualise not only attractors but also their basins.
In the realm of asynchronous, non-deterministic modeling not only is the repertoire of software even more limited, but also the formal notions for basins of attraction are often lacking.
In this setting, the difficulty both for theory and computation arises from the fact that states may be elements of several distinct basins.
In this paper we address this topic by partitioning the state space into sets that are committed to the same attractors.
These commitment sets can easily be generalised to sets that are equivalent w.r.t. the long-term behaviours of pre-selected nodes which leads us to the notions of markers and phenotypes which we illustrate in a case study on bladder tumorigenesis.
For every concept we propose equivalent CTL model checking queries and an extension of the state of the art model checking software NuSMV is made available that is capable of computing the respective sets.
All notions are fully integrated as three new modules in our Python package PyBoolNet, including functions for visualising the basins, commitment sets and phenotypes as quotient graphs and pie charts.









§ BACKGROUND



Boolean networks offer an intuitive approach to simulating the dynamics of interaction networks.
In cell biology these are usually gene regulatory or signal transduction networks.
Each biological component is modelled by a Boolean variable that can only switch between two states, true and false.
Depending on the component, these states represent whether a protein is present at high or low concentration levels,
a gene is being transcribed at or above its base rate, a signalling molecule is phosphorylated or not, and so on.
The conditions under which variables switch on and off are specified by Boolean functions that involve other variables of the network.
Together with a transition relation, which determines how many and which components are allowed to change during a state transition,
Boolean networks can be used to predict transient and long term activity profiles of the involved components.
Although the restriction to discrete states is arguably an over-simplification for many processes,
Boolean networks have frequently been shown to be fruitful to predict drug targets, phenotypes, novel interactions, etc.
This paper is particularly relevant to models whose state transitions are non-deterministic.
A prominent example are the so-called Thomas networks (after René Thomas, see e.g. <cit.>) where changes in state variables are modelled as processes that are not synchronised.
Asynchronous transition systems are especially applicable when details about the precise mechanics of the processes, e.g. the activation and inhibition delays, are unknown or uncertain.
The results of this paper are particularly relevant to non-deterministic transition relations but also work in the deterministic case.
For a detailed introduction to Boolean networks and qualitative modelling more generally see <cit.>.

Attractor detection is a frequent starting point for testing the consistency of models with data and also for obtaining model predictions, see e.g. <cit.>.
An extension of this approach considers not only the location of attractors but also their basins of attraction.
These regions and transitions between them determine the likelihood of reaching particular attractors and not others.
From a practical, wet-lab point of view, the motivation for studying basins and attractor commitment sets is to find out how the underlying decision processes are implemented within the interaction networks of living cells.
Computational models and the formalisation of notions like basin, commitment set and phenotype add to this research by making hypothesis testable in silico.
When these models are validated they may then be used to measure the robustness of interaction networks with respect to knock-outs or knock-ins and other perturbations.
Due to the high computational cost of analysing the complete state transition graph for larger networks, such studies often use aggregated results from simulations  <cit.>. 
While analytic results are mostly focused on Boolean networks with deterministic updates <cit.>,
there are still some algorithms for computing the basins of Boolean networks with asynchronous updates <cit.>.
These all have to deal with the inherent difficulties in analysing asynchronous state transition graphs, which in turn results in restricted scalability.

In general, the computation of basins of attraction is challenging because the size of the state space grows exponentially with the number of variables of the model.
A competitive approach for dealing with large transition systems is to use so-called symbolic model checking algorithms <cit.>.
Model checking is a method from computer science to decide whether formal hardware designs satisfy desired specifications.
The designs are formalised as so-called Kripke structures, which are labelled, directed graphs,
and the specifications to be verified are translated into statements of a chosen temporal logic.
Computation Tree Logic (CTL) is an example of such a logic.
It allows statements about the alternative futures originating in states <cit.>.
These are particularly useful when dealing with properties related to attractors.
NuSMV is an open source symbolic model checking tool allowing for the efficient analysis of transition systems <cit.>.

It has been observed that the state transition graphs (STGs) of Boolean networks are Kripke structures and that model checking is therefore applicable to Boolean networks <cit.>.
There are custom model checking tools that are specifically developed for biological networks, e.g. <cit.>,
but their additional capabilities may come at the cost of efficiency in more basic queries.
Still, we believe that for the problem of computing basins of attractors and related objects state of the art software designed to answer reachability queries exactly such as the model checking software NuSMV provides the ideal basis for implementation.

In this paper we, on the one hand, formally develop several notions related to basins of attraction and, on the other hand, provide efficient tools to compute these objects in practice using a model checking approach.
Namely, we created NuSMV-a, an extension of NuSMV 2.6, which implements model checking with accepting states.
We complement the theoretical part in later sections by giving an impression of the tool efficiency and scalability and by considering a case study on bladder tumorigenesis,
focusing on analysis of different phenotypes of the cancer cells.

We aim at both clear and intuitive descriptions, further technical details can be found in the supplementary material.
In addition, scripts and tutorials to reproduce the results are available, integrated in our PyBoolNet software <cit.>.




§ MODEL CHECKING WITH ACCEPTING STATES

CTL model checking was originally designed to answer queries by either yes or no depending on whether the transition system satisfies the query or not.
When model checking is used commercially, i.e., for the verification of software and hardware designs,
it turned out to be useful to analyse the reasons for why a design does not satisfy a specification.
This lead to the development of algorithms that generate counterexamples when queries are false <cit.>.
When model checking is applied to systems biology a similar analysis is desirable.
In particular, the theory which is developed in the subsequent sections requires the capability of computing the exact set of states that satisfies a query, an idea that is more basic then computing counterexamples.
The accepting states of a CTL formula φ and STG (S,→) are defined by

    𝐴𝑐𝑐𝑒𝑝𝑡(S,→,φ) := {x∈ S | xφ}

where xφ means that φ holds in state x of the transition system (S,→).
For more details on model checking and CTL, see the supplementary material.
Note that the accepting states are defined without initial states.
By definition, see <cit.>, a transition system (S,→) satisfies a model checking query with CTL formula φ
and initial states I⊆ S iff I⊆𝐴𝑐𝑐𝑒𝑝𝑡(S,→,φ).

For our purposes, the accepting states of one query are often re-used in the CTL formula of a subsequent query.
We therefore need to be able to encode sets of states as Boolean expressions.
Since the state space is already defined in terms of states of Boolean variables,
we may use a disjunctive normal form (DNF) in which each term references a single state of the set.
Clearly, this is not a very efficient encoding and there are usually other, shorter expressions that also represent the set.
We denote by 𝐸𝑛𝑐(X) an arbitrary Boolean expression that represents X⊆ S.
Therefore, 𝐸𝑛𝑐(𝐴𝑐𝑐𝑒𝑝𝑡(S,→,φ)) is an expression that represents the accepting states of φ in (S,→).





§ BASINS OF ATTRACTION

Basins of attraction describe regions of state space that lead into specific attractors.
Large basins suggest a certain robustness of their attractors against perturbations, see e.g. <cit.>, while small basins suggest unlikely attractors.
The study of basins is therefore essential for understanding the likely long-term behaviours of Boolean networks.

The basins of deterministic and non-deterministic STGs are very different
because states in deterministic STGs can always reach exactly one attractor while states in non-deterministic STGs can usually reach several attractors.
It has therefore been suggested, for example in <cit.>, to distinguish between weak and strong basins of attraction
where the weak basin contains all states that can reach a given attractor while the strong basin contains all states that can only reach a given attractor.
For an attractor A⊆ S and an arbitrary state a∈ A, the weak basin of an attraction is defined by

    𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(S,→,a) := {x∈ S |∃π∈𝑃𝑎𝑡ℎ𝑠(x): a∈π},

where 𝑃𝑎𝑡ℎ𝑠(x) is the set of all paths starting from x in the (S,→) and a∈π means that a appears somewhere along π.
Note that, as attractors are strongly connected, it is irrelevant which a∈ A is chosen.
Biologically, these basins represent the weakest form of commitment to an attractor and the system might never actually reach the given attractor,
either because there are other attractors that are more likely or because the system remains in an intermediate equilibria.

To define the strong basin one may be tempted to require that all paths lead to the attractor.
But this definition is stronger than the notion that a state can only reach A.
It excludes states that are on cycles, even if they can, in terms of attractors, only reach A.
Nevertheless this type of basin is interesting and since it is cycle-free, we call it 𝐶𝑦𝑐𝐹𝑟𝑒𝑒𝐵𝑎𝑠𝑖𝑛(S,→,A) and define it by

    {x∈ S |∀π∈𝑃𝑎𝑡ℎ𝑠(x): ∃ a∈ A: a∈π}.

Note that this definition is in terms of A instead of an arbitrary a∈ A.
The reason is that in specific cases it makes a difference which a∈ A is chosen, e.g., if A is a complex attractor that contains a cycle that does not visit all states of A.
These basins form the strictest notion of attraction, the inevitable fall into the given attractor without the possibility of intermediate oscillations or the potential of reaching another attractor.

The strong basin in the usual sense of all states that can only reach A is defined by called 𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛(S,→,a) and defined by

    {x∈ S |∀π∈𝑃𝑎𝑡ℎ𝑠(x): ∀ y∈π: ∃π'∈𝑃𝑎𝑡ℎ𝑠(y):  a∈π'}.

As the system transitions from a weak to a strong basin it looses the potential of reaching certain other attractors and hence these transitions hint at hidden decision processes within the changing activity patterns.
These three notions are illustrated for a toy transition system in Fig.<ref>.
The Boolean equations for the toy transition system were created using iBoolNet <cit.>.



To apply symbolic model checking to the computation of weak and strong basins we need CTL queries that capture these definitions.
We only need three operators, namely  the exists finally operator 𝐄𝐅, the always finally operator 𝐀𝐅 and the all globally operator 𝐀𝐆.
These operators can be linked to the path-based definitions of the basins, since a state x satisfies 𝐄𝐅(M) for a state set M if there exists a path starting in x leading to some state in M,
similarly 𝐀𝐅(M) encodes a condition for all paths starting in x and the combination 𝐀𝐆(𝐄𝐅(M)) is satisfied if from every state reachable from x there exists a path to a state in M.
This last condition would in particular not be true for a state x that can reach an attractor M and in addition a second attractor M', since then the states of M' are reachable from x, but M is not reachable from M'.
The formal definitions can be found in the supplementary material. These considerations lead us to the following CTL queries

 𝐴𝑐𝑐𝑒𝑝𝑡(S,→,α),   where  α:=𝐄𝐅(𝐸𝑛𝑐(a))

 𝐴𝑐𝑐𝑒𝑝𝑡(S,→,β),   where  β:=𝐀𝐅(𝐸𝑛𝑐(A))

 𝐴𝑐𝑐𝑒𝑝𝑡(S,→,γ),   where  γ:=𝐀𝐆(𝐄𝐅(𝐸𝑛𝑐(a)))

for 𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(a), 𝐶𝑦𝑐𝐹𝑟𝑒𝑒𝐵𝑎𝑠𝑖𝑛(A) and 𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛(a), respectively.
Note that the three types of basin satisfy the inclusions

    A ⊆𝐶𝑦𝑐𝐹𝑟𝑒𝑒𝐵𝑎𝑠𝑖𝑛(A) ⊆𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛(a) ⊆𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(a)



To visualise the relative sizes of the basins we suggest to use pie charts and stacked bar plots.
The bar plots highlight the inclusion relationship while the pie charts give an impression of the relative sizes of the basins.
These three notions are illustrated for the toy transition system in Fig.<ref>.



The notion of basin can be generalised from single to multiple attractors simply by adjusting the state formulas accordingly.
Suppose we are interested in the basins of k attractors A_1,…,A_k with representative states a_i∈ A_i.
With X:={a_1,…,a_k} and Y:=A_1∪…∪ A_k the CTL queries for the generalised basins are obtained by using X and Y instead of a and A in the definitions above:
𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(X), 𝐶𝑦𝑐𝐹𝑟𝑒𝑒𝐵𝑎𝑠𝑖𝑛(Y) and 𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛(X).
Note that while the distributive law 𝐄𝐅(ϕ_1∨ϕ_2)=𝐄𝐅(ϕ_1)∨𝐄𝐅(ϕ_2) for existential CTL operators implies
that the weak basin of several attractors is the union of the individual weak basins:

    𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛({a_i,a_j}) = 𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(a_i) ∪𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(a_j)

the same does not hold for the strong and cycle-free basins.
Instead, the following inclusions hold:

 𝐶𝑦𝑐𝐹𝑟𝑒𝑒𝐵𝑎𝑠𝑖𝑛(A_i ∪A_j)    ⊇𝐶𝑦𝑐𝐹𝑟𝑒𝑒𝐵𝑎𝑠𝑖𝑛(A_i)

 𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛({a_i,a_j})      ⊇𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛({a_i})

The reason for the inclusion instead of equality is that a strong basin of two attractors may contain a state that can reach both attractors.
Such a state can not belong to the individual strong basins.
The same example proves that the inclusion may be strict for cycle-free basins.

Note that all definitions above work even if not all attractors are known.
In practice we may therefore compute, for example, the strong basin of a steady state even if we do not know whether there are additional attractors.




§ COMMITMENT SETS

A different approach to understanding attractors and regions in state space that lead to them is via commitment sets.
Commitment sets partition the state space into classes that can reach the exact same attractors.
For ease of referencing specific commitment sets, let (A_1,…,A_k) be a sequence of k attractors with representatives (a_1,…,a_k)
and let I⊆{1,…,k} be a non-empty subset of indices.
The states that are committed to the attractors {A_i | i∈ I} are defined by

    𝐶𝑜𝑚𝑆𝑒𝑡(I) := {x∈ S | (∃π∈𝑃𝑎𝑡ℎ𝑠(x):  a_i∈π) ⇔ (i∈ I)}.

Commitment sets are the finest resolution of state space with respect to the potential of reaching certain attractors.
An analysis of the similarities of states that belong to the same set is biologically as relevant as understanding the differences between states of different commitment sets.
In particular, it remains an open question to understand the decision processes that are realised by transitions between commitment sets, see also <cit.>.

To encode this in CTL we observe that a state is committed to the attractors indexed by I if and only if there is a path from x to a_i for every i∈ I
and x is in the strong basin of {A_i | i∈ I}, i.e., can not reach other attractors.
A possible CTL query for 𝐶𝑜𝑚𝑆𝑒𝑡(I) is therefore 𝐴𝑐𝑐𝑒𝑝𝑡(S,→,δ), where


    δ:=⋀_i∈ I𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(a_i) ∧𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛({a_j | j∈ I}).


Note that, as before, commitment sets can be computed without knowing all attractors of the network.
Furthermore, the same queries can be used to compute the basins and commitment sets of subspaces and in particular minimal trap spaces, see <cit.>,
which are often in a one-to-one correspondence with the attractors of a network, see <cit.>.

For small networks with up to 15 nodes or so we may visualise commitment sets by different node colours.
An example is given in Fig.<ref>.
For larger networks we suggest to visualise the commitment sets and transitions between them by a "graph of commitment sets" in which the nodes are commitment sets and there are edges between the nodes if there are transitions between the corresponding commitment sets in the underlying STG.
More precisely, we construct the quotient graph of the STG, induced by the commitment sets.
Recall from basic graph theory that the quotient graph of a digraph G=(N,→), induced by a partition N_1∪…∪ N_m=N of its nodes N,
is the digraph Q whose nodes are the blocks N_i and there is an edge between
two blocks N_i→ N_j if there are nodes n_1∈ N_i and n_2∈ N_j with and edge n_1→ n_2 in G.

In the context of commitment sets we call the quotient graph of an STG induced by the commitment sets its commitment diagram.
An example is given in Fig.<ref>.

As for the nodes, we compute the edges of the commitment diagram by answering CTL model checking queries, using the "exists next" operator (𝐄𝐗) which is satisfied for a state x if the given formula holds in a state that is reachable from x by a single transition.
Given two commitment sets 𝐶𝑜𝑚𝑆𝑒𝑡(I) and 𝐶𝑜𝑚𝑆𝑒𝑡(J) with I,J⊂{1,…,k}, there is an edge 𝐶𝑜𝑚𝑆𝑒𝑡(I)→𝐶𝑜𝑚𝑆𝑒𝑡(J) if
there are states x_1∈𝐶𝑜𝑚𝑆𝑒𝑡(I) and x_2∈𝐶𝑜𝑚𝑆𝑒𝑡(J) such that x_1→ x_2 in the STG.
Using model checking we first compute


    𝐴𝑐𝑐𝑒𝑝𝑡(S,→,ζ),   where  ζ:=𝐄𝐗(𝐸𝑛𝑐(𝐶𝑜𝑚𝑆𝑒𝑡(J))).


and note that the edge 𝐶𝑜𝑚𝑆𝑒𝑡(I)→𝐶𝑜𝑚𝑆𝑒𝑡(J) exists if

    𝐶𝑜𝑚𝑆𝑒𝑡(I) ∩𝐴𝑐𝑐𝑒𝑝𝑡(S,→,ζ) ≠∅.

Note that an edge can only exist if I⊋ J because the set of reachable attractors is decreasing along transitions in an STG.













§ PHENOTYPES

In this section we extend the notion of commitment sets to make it more suited to application.
In experimental studies rarely all components of a system can be measured at the same time, if even at all.
So when using attractors for model validation and prediction an obvious drawback is that the experimental data yields only partial descriptions that might be valid for several attractors.
To be more faithful to this scenario, attractors can be grouped together according to the values of so-called marker components.
The resulting object is called a phenotype.
More specifically, a phenotype is a pattern that specifies the long-term behaviours of a set of given markers with respect to three possible outcomes: inhibited=0, expressed=1 and oscillating=*.
A phenotype analysis can help to understand fundamental system properties by simplifying the attractor landscape of the system.
Its main advantage in application, however, is that the measurable components naturally form a set of markers that link the predictions of the Boolean network with observed expression data <cit.>.

We formalise phenotypes as three-valued functions p:U→{0,1,*} from the set of markers U⊆ V to the set of possible outcomes {0,1,*}.
The phenotype p_A of an attractor A⊆ S is defined by

    p_A(u)=0:   ∀ a∈ A:  a(u)=0 
     1:   ∀ a∈ A:  a(u)=1 
    :   else

Given a set of markers U we say that a phenotype p exists if there is an attractor A⊆ S such that p_A=p.
Note that there may be several attractors with identical phenotypes, p_A=p_A', but every attractor has a unique phenotype (for fixed markers).

Analogous to the commitment diagram we define the phenotype diagram by partitioning the state space into subsets that are committed to identical phenotypes.
Let p_1,…,p_k be the phenotypes of a given STG and markers.
For each phenotype p_i, let A_i_1,…,A_i_l be all corresponding attractors with representative states a_i_j∈ A_i_j and denote by P_i:={a_i_1,…,a_i_l} the union of the representative states.
We call the commitment set for the phenotypes specified by I⊆{1,…,k} the phenotype set.
It is called 𝑃ℎ𝑒𝑛𝑜𝑆𝑒𝑡(I) and defined by

    {x∈ S | (∃π∈𝑃𝑎𝑡ℎ𝑠(x): ∃ a∈ P_i:  a∈π) ⇔ (i∈ I)}.

It can be computed using model checking with accepting states 𝐴𝑐𝑐𝑒𝑝𝑡(S,→,η) where


    η:=⋀_i∈ I𝑊𝑒𝑎𝑘𝐵𝑎𝑠𝑖𝑛(P_i) ∧𝑆𝑡𝑟𝑜𝑛𝑔𝐵𝑎𝑠𝑖𝑛(⋃_i∈ I P_i).

Note the similarity with the query for commitment sets.
The only difference is that we use sets of representatives P_i instead of individual representatives a_i.
As before, the existence of edges 𝑃ℎ𝑒𝑛𝑜𝑆𝑒𝑡(I)→𝑃ℎ𝑒𝑛𝑜𝑆𝑒𝑡(J) between phenotype sets is determined via the quotient graph relation.




§ NUSMV WITH ACCEPTING STATES

To compute the basin, commitment and phenotype diagrams, we developed an extension of NuSMV that is capable of returning the accepting states of a transition system and CTL formula.
It is based on NuSMV 2.6, called NuSMV-a ("a" for "accepting states"), available from <cit.> and fully integrated in PyBoolNet starting from version 2.0 <cit.>.
In order to decide whether a CTL formula φ is true or false,
NuSMV already computes the accepting states in the form of a BDD that represents 𝐴𝑐𝑐𝑒𝑝𝑡(S,→,φ).
We only had to make this information accessible and implement a suitable format in which to return the accepting states.
We decided to use an output format called factored form that is implemented in the function  (CUDD 3.1.4).
Although the resulting Boolean expressions become quickly unreadable for humans they are in NuSMV syntax and can therefore be used as sub-formulas of subsequent CTL queries.
To obtain minimised expressions that may be easier for humans to read we extended PyBoolNet with a module called  and
specifically the function  that uses the program espresso for minimising Boolean expressions <cit.>.

As of version 2.0, PyBoolNet comes with NuSMV-a and the function

    
 of the package  returns
a Python dictionary that contains all formulas and cardinalities as strings and integers, respectively.
The PyBoolNet modules for computing the various basins, the commitment sets and the phenotypes all require this new model checking function.
For tutorials and examples consult the PyBoolNet manual and project homepage.




§ EFFICIENCY AND SCALABILITY

This section presents a table of running times for computing the accepting states of four Boolean networks of increasing size:
 <cit.> with n=12,
 <cit.> with n=24,
 <cit.> with n=53 and
 <cit.> with n=103.
The networks are available in the PyBoolNet model repository.
We generated three CTL queries of increasing complexity for each network: a weak basin query, a strong basin query and a commitment set query.
For state formulas we compute two attractor states a_i∈ A_i for i=1,2 for each network using minimal trap spaces and used a conjunctive encoding.
We compared the running times of NuSMV-a with Antelope (Analysis of Networks through Temporal-Logic Specifications) the only other software we found for computing accepting states of CTL queries.
Antelope is a model checking software specifically designed for biological interaction networks, see <cit.>.
Its model checking algorithms are based on "hybrid CTL", an extension of standard CTL that allows, for example, the specification of cycles of a given length.
We used a 64-bit Ubuntu 16.04 LTS desktop pc with 7.7 GiB memory and 4 CPUs with 2.7 GHz and the NuSMV-a command .
The results in Tab. <ref> show that NuSMV can handle networks with up to around 50 variables in reasonable times in the order of minutes.
Antelope can answer queries in reasonable time for about 20 variables.












§ CASE STUDY: BLADDER TUMORIGENESIS

In this case study we investigate the bladder tumorigenesis model of <cit.> and compare the results the trajectory-based studies that are also available in <cit.>.
The tumorigenesis model has four inputs which represent the cellular context of the network in terms of
DNA damage, the growth factors EGFR and FGFR3 and growth inhibitors. Three output components, Growth
Arrest, Proliferation and Apoptosis determine the response of the cell to a given context.
The model was constructed using the software GINsim <cit.>.
There are five multi-valued components: ATM, Apoptosis, CHEK1_2, E2F1 and E2F3 each with the activity levels 0=low, 1=medium and 2=high.
Since we use NuSMV-a as a plugin of PyBoolNet which is designed for Boolean networks we worked with the equivalent Booleanised version in which multi-valued components x are replaced by two Boolean components x_medium and x_high.
The update functions are obtained by the so-called van Ham mapping which has been shown to result in an equivalent transition system <cit.>.
PyBoolNet detects variables that represent the same multi-valued variable and removes inadmissable combinations of values from the initial states of model checking queries.
Meaningless states, e.g., x_medium=0 and x_high=1, that are artefacts of the Booleanisation do therefore not contribute to the basins, commitment sets and phenotype sets.
The Booleanisation was done using GINsim and the Boolean model is available from the PyBoolNet model repository under the name remy_rumorigenesis.

The analysis of the wild-type multi-valued model consists of a list of the observed attractors and their phenotypes, see Table S2 of the supplementary materials of <cit.>, here referred to as Table S2,
and pie charts that visualise the reachability of each phenotype under different conditions, e.g. the wildtype but also various knock-out mutants.
The reachabilities of phenotypes are defined as the percentage of simulations (random walks on the asynchronous state transition graph) that exhibit a certain phenotype.
The authors ran 10,000 simulations for each of the 16 input configurations where the remaining variables are randomly initialised.

We first confirmed that the attractors, 20 steady states and 5 cyclic attractors, of Table S2 are identical to the ones we find using the minimal trap space method <cit.> for our Boolean model.
The 25 attractors exhibit four phenotypes for the markers Proliferation, Apoptosis_medium, Apoptosis_high and Growth_arrest.
The phenotypes are listed in Tab. <ref>.

The basins of attraction are visualised as a bar plot in Fig. <ref>.
Since the model has four inputs, its state space is divided into 16 disconnected regions each making up 1/16=6.25% of all states.
The maximal height of a bar is therefore at 6.25.
It can be deduced that the first 9 attractors are the only ones in the state space region defined by their respective input combination because the whole region belongs to their strong basins.
The remaining 16 attractors share their input combination with another attractor and the majority of states belong to more than one basin.



The commitment diagram also reflects the division into 16 disconnected regions.
It is divided into 16 disconnected diagrams of three different shapes where the shape depends only on how many attractors there are in the respective region.
An example of each is given in Fig. <ref>.
If an input combination contains more than one attractor, Fig. <ref> (a) and (b), then the majority of states, about 75%, can reach all attractors.



We then compared the phenotype reachability analysis of <cit.> with our notion of phenotype commitment.
The reachability analysis computes the percentage of simulations that exhibit a given phenotype.
The results are visualised as a pie chart with one slice per phenotype, see Fig. <ref> (a) (taken from <cit.>).
Our state-based analysis, which uses model checking with accepting states, finds that there are five different phenotype commitment sets:
one for each phenotype but also one additional slice that represents states in the basin of several phenotypes, namely    and , see Fig. <ref> (b).
To make the comparison between the trajectory-based analysis of <cit.> and our state-based analysis easier we also visualise the commitment sets using a pie chart and use the same colours.
Looking at both charts suggests that the majority of the states not yet committed to a single attractor (namely 10/12%) lead to  but neither chart alone holds the full information of the probability of reaching the phenotypes and the sizes of the commitment sets.
The two types of analysis can be seen as complementing each other.

















§ DISCUSSION


This paper puts the notions of a basin of attraction for non-deterministic transition systems on firm mathematical ground.
We propose definitions of varying strictness and focus that allow to partition state space in a way providing insight in the sequences of both reversible and non-reversible decisions ultimately leading the system into an attractor.
For each introduced set and its corresponding reachability properties, we developed means of visualisation, either in the form of a diagram or a pie chart.

On the computational side we define what CTL model checking with accepting states is and offer an extension to NuSMV that is capable of returning the accepting states of a query.
This extension, NuSMV-a, as well as the computation and visualisation tools for all notions introduced here are now part of PyBoolNet which offers a high level interface to model checking Boolean networks with many functions for model and query generation already implemented.
The comparison of NuSMV-a with Antelope suggests that computing accepting states with NuSMV-a is competitive and networks with more than 50 nodes can be handled on a scale of minutes.

Using the case study, we illustrated the various diagrams for a molecular network and also to confirm the correctness of the results in comparison with the study of Remy et al., e.g., concerning the calculated attractors and phenotypes.
We showed that our approach yields a useful complementary view on assessing the likelihood of reaching particular attractors and the corresponding variability inherent in the system.

We believe based on our results that in the context of non-deterministic transition systems a range of notions concerning attractor reachability are of interest.
The cycle-free basins springing from the natural condition that all paths originating in a state should lead to the same attractor show that further distinctions concerning the asymptotic behavior of trajectories are possible and useful.
The question in how far non-terminal cycles in the state transition graph can play a similar role to attractors lead to questions about the stability of such cycles.
Progress concerning this issue in turn could potentially be highly meaningful for application, since it may allow to uncover robustly maintainable system behaviors in the models that are overlooked when only focusing on the attractors.
A possible way to address this question mathematically and computationally may be in combining our state based analysis of commitment sets with trajectory-based analysis as we have outlined in the case study for the reachability analysis of phenotypes.

In a similar step to higher resolution, the commitment sets can be further refined by considering not only which attractors are reachable but also which other commitment sets, i.e.,
two states may both be committed to the same attractors but they traverse different commitment sets on their way to the attractors.
It is certainly of mathematical interest to find the rule allowing to link the properties of the resulting graphs on the different resolution levels, which in turn links to the most exciting open question that results from this research is the question of characterising the control points and decision processes that occur at the border between commitment sets.
Which changes in activity of which variables are responsible for the loss of a previously attainable attractor?
What are minimal interventions that restore the reachability?
Questions along these lines are crucial for understanding the dynamics of qualitative models and in practice for controlling cell differentiation and cell reprogramming.




§ ACKNOWLEDGMENT



We thank the participants of the discussion about basins of attraction at Luminy University in Marseille in 2016:
Elisabeth Rémy, Claudine Chaouiya, Laurence Calzone, Denis Thieffry, Aurélien Naldi, Pedro Monteiro and Laurent Tournier.
The work was partially funded by the German Federal Ministry of Education and Research (BMBF), grant no. 0316195.



unsrt


