A quick tutorial showing the features of BoNesis for the synthesis of Most Permissive Boolean Networks from network architecture and dynamical constraints.
import bonesis
import pandas as pd
from colomoto_jupyter import tabulate
The synthesis uses two inputs:
Let us define an influence graph from a list of pairwise interactions, with a sign.
influences = [
("Pax6","Pax6",dict(sign=1)),
("Pax6","Hes5",dict(sign=1)),
("Pax6","Mash1",dict(sign=1)),
("Hes5","Mash1",dict(sign=-1)),
("Hes5","Scl",dict(sign=1)),
("Hes5","Olig2",dict(sign=1)),
("Hes5","Stat3",dict(sign=1)),
("Mash1","Hes5",dict(sign=-1)),
("Mash1","Zic1",dict(sign=1)),
("Mash1","Brn2",dict(sign=1)),
("Zic1","Tuj1",dict(sign=1)),
("Brn2","Tuj1",dict(sign=1)),
("Scl","Olig2",dict(sign=-1)),
("Scl","Stat3",dict(sign=1)),
("Olig2","Scl",dict(sign=-1)),
("Olig2","Myt1L",dict(sign=1)),
("Olig2","Sox8",dict(sign=1)),
("Olig2","Brn2",dict(sign=-1)),
("Stat3","Aldh1L1",dict(sign=1)),
("Myt1L","Tuj1",dict(sign=1)),
]
dom1 = bonesis.InfluenceGraph(influences)
dom1
# computing graph layout...
Here, dom1
delimits any BN that uses at most the listed influences, with the right sign. Thus, some solutions may use only a subset of this influence graph.
If you want to enforce BNs using all the given influences, use the option exact=True
:
dom2 = bonesis.InfluenceGraph(influences, exact=True)
For influence graphs with large in-degrees, it is necessary to specify a bound on the number of clauses in the disjunction normal form (DNF) of the BNs with the maxclause
argument. See help(bonesis.InfluenceGraph)
for other options.
They are specified by a Python dictionnary associating observation names to observed states of a subset of nodes:
data = {
"zero": {n: 0 for n in dom1}, # all nodes are 0
"init": {n: 1 if n == "Pax6" else 0 for n in dom1}, # all nodes are 0 but Pax6
"tM": {"Pax6": 1, "Tuj1": 0, "Scl": 0, "Aldh1L1": 0, "Olig2": 0, "Sox8": 0},
"fT": {"Pax6": 1, "Tuj1": 1, "Brn2": 1, "Zic1": 1, "Aldh1L1": 0, "Sox8": 0},
"tO": {"Pax6": 1, "Tuj1": 0 ,"Scl": 0, "Aldh1L1": 0, "Olig2": 1, "Sox8": 0},
"fMS": {"Pax6": 1, "Tuj1": 0, "Zic1": 0, "Brn2": 0, "Aldh1L1": 0, "Sox8": 1},
"tS": {"Pax6": 1, "Tuj1": 0, "Scl": 1, "Aldh1L1": 0, "Olig2": 0, "Sox8": 0},
"fA": {"Pax6": 1, "Tuj1": 0, "Zic1": 0, "Brn2": 0, "Aldh1L1": 1, "Sox8": 0},
}
pd.DataFrame.from_dict(data, orient="index").fillna('')
Pax6 | Hes5 | Mash1 | Scl | Olig2 | Stat3 | Zic1 | Brn2 | Tuj1 | Myt1L | Sox8 | Aldh1L1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
zero | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
init | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
tM | 1 | 0 | 0 | 0 | 0 | 0 | ||||||
fT | 1 | 1 | 1 | 1 | 0 | 0 | ||||||
tO | 1 | 0 | 1 | 0 | 0 | 0 | ||||||
fMS | 1 | 0 | 0 | 0 | 1 | 0 | ||||||
tS | 1 | 1 | 0 | 0 | 0 | 0 | ||||||
fA | 1 | 0 | 0 | 0 | 0 | 1 |
bo = bonesis.BoNesis(dom1, data)
The data
dictionnary specifies observations that can be used to constraint configurations (or states) of the network.
There are two shortcuts for binding a configuration to an observation:
~bo.obs("A")
is a unique pre-defined configuration bound to the observation "A"
;
+bo.obs("A")
returns a new configuration bound to "A"
. Thus in the following code:
cfg1 = ~bo.obs("A")
cfg2 = ~bo.obs("A")
cfg3 = +bo.obs("A")
cfg4 = +bo.obs("A")
cfg1
and cfg2
refers to the same configuration; whereas cfg3
and cfg4
may be different.
We detail two kind of attractor constraints: fixed points and trap spaces. Both are specified with the fixed
predicate, which, depending on the argument will enforce the existence of one of the two kinds of attractor.
When giving a configuration as argument, fixed
ensures that the configuration is a fixed point:
bo.fixed(~bo.obs("fA"))
bo.fixed(~bo.obs("fMS"));
A trap space specification is given by an observation, which enforces that all the nodes in the given observations can never change of value (thus any reachable attractor have these nodes fixed):
fT_tp = bo.fixed(bo.obs("fT"))
~bo.obs("init") >= ~bo.obs("tM") >= fT_tp
~bo.obs("init") >= ~bo.obs("tO") >= ~bo.obs("fMS")
~bo.obs("init") >= ~bo.obs("tS") >= ~bo.obs("fA");
They can be specified using the nonreach
function, or equivalently the /
operator between two configurations. The right-hand side can also be a fixed
constraint.
~bo.obs("zero") / fT_tp
~bo.obs("zero") / ~bo.obs("fMS")
~bo.obs("zero") / ~bo.obs("fA");
Enumerations of solutions are done through iterators. The basic one being the boolean_networks
which returns mpbn.MPBooleanNetwork
objects.
for bn in bo.boolean_networks(limit=2): # limit is optional
print(bn)
Grounding...done in 0.0s Aldh1L1 <- Stat3 Brn2 <- Mash1 Hes5 <- !Mash1&Pax6 Mash1 <- !Hes5&Pax6 Myt1L <- 0 Olig2 <- Hes5&!Scl Pax6 <- Pax6 Scl <- Hes5&!Olig2 Sox8 <- Olig2 Stat3 <- Hes5&Scl Tuj1 <- Brn2 Zic1 <- Mash1 Aldh1L1 <- Stat3 Brn2 <- Mash1 Hes5 <- !Mash1 Mash1 <- !Hes5&Pax6 Myt1L <- 1 Olig2 <- Hes5&!Scl Pax6 <- Pax6 Scl <- Hes5&!Olig2 Sox8 <- Olig2 Stat3 <- Hes5&Scl Tuj1 <- Zic1 Zic1 <- Mash1
Display as a table:
solutions = list(bo.boolean_networks())
pd.DataFrame(solutions)
Grounding...done in 0.0s
Aldh1L1 | Brn2 | Hes5 | Mash1 | Myt1L | Olig2 | Pax6 | Scl | Sox8 | Stat3 | Tuj1 | Zic1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Stat3 | Mash1 | !Mash1 | !Hes5&Pax6 | 1 | Hes5&!Scl | Pax6 | Hes5&!Olig2 | Olig2 | Scl | Brn2&Myt1L | Mash1 |
1 | Stat3 | Mash1 | !Mash1 | !Hes5&Pax6 | 0 | Hes5&!Scl | Pax6 | Hes5&!Olig2 | Olig2 | Scl | Zic1 | Mash1 |
2 | Stat3 | Mash1 | !Mash1 | !Hes5 | 1 | Hes5&!Scl | Pax6 | Hes5&!Olig2 | Olig2 | Scl | Brn2&Zic1 | Mash1 |
3 | Stat3 | Mash1 | !Mash1&Pax6 | !Hes5&Pax6 | Olig2 | Hes5&!Scl | Pax6 | Hes5&!Olig2 | Olig2 | Scl | Brn2&Zic1 | Mash1 |
4 | Stat3 | Mash1 | !Mash1&Pax6 | !Hes5&Pax6 | Olig2 | Hes5&!Scl | Pax6 | Hes5&!Olig2 | Olig2 | Scl | (Brn2&Myt1L)|(Brn2&Zic1) | Mash1 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
1115 | Stat3 | Mash1&!Olig2 | !Mash1&Pax6 | !Hes5 | 0 | Hes5&!Scl | Pax6 | !Olig2 | Olig2 | Hes5&Scl | Brn2|Zic1 | Mash1 |
1116 | Stat3 | Mash1 | !Mash1&Pax6 | !Hes5 | 0 | Hes5&!Scl | Pax6 | !Olig2 | Olig2 | Hes5&Scl | Brn2&Zic1 | Mash1 |
1117 | Stat3 | Mash1&!Olig2 | !Mash1&Pax6 | !Hes5 | 0 | Hes5&!Scl | Pax6 | !Olig2 | Olig2 | Hes5&Scl | Brn2&Zic1 | Mash1 |
1118 | Stat3 | Mash1 | !Mash1 | !Hes5&Pax6 | 0 | Hes5&!Scl | Pax6 | !Olig2 | Olig2 | Hes5&Scl | Brn2|Zic1 | Mash1 |
1119 | Stat3 | Mash1&!Olig2 | !Mash1 | !Hes5&Pax6 | 0 | Hes5&!Scl | Pax6 | !Olig2 | Olig2 | Hes5&Scl | Brn2|Zic1 | Mash1 |
1120 rows × 12 columns
len(solutions)
1120
Resulting objects are MPBooleanNetwor
from the mpbn
Python library. One can then compute reachability and attractor properties directly:
pd.DataFrame(solutions[0].attractors())
Aldh1L1 | Brn2 | Hes5 | Mash1 | Myt1L | Olig2 | Pax6 | Scl | Sox8 | Stat3 | Tuj1 | Zic1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 0 | 0 |
1 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 |
2 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 |
3 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 |
4 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 1 |
pd.DataFrame(solutions[0].attractors(reachable_from=data["init"]))
Zic1 | Tuj1 | Stat3 | Sox8 | Scl | Pax6 | Olig2 | Myt1L | Mash1 | Hes5 | Brn2 | Aldh1L1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 |
1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 0 |
2 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 1 |
bo.all_fixpoints({bo.obs(obs) for obs in ["fA", "fMS", "fT", "zero"]});
The following constraint ensures that any fixed point reachable from the configuration on left-hand side has to match with at least one of the given observation.
~bo.obs("init") >> "fixpoints" ^ {bo.obs(obs) for obs in ["fA", "fMS", "fT"]};
bo.boolean_networks().count()
Grounding...done in 0.0s
88
To better understand the composition of the different solutions, one can project the solutions on each node: given a node A, it enumerates all the Boolean functions for A that are used in at least one full solution.
The projected solutions can be accessed from the following object:
projs = bo.local_functions()
Grounding...done in 0.0s
The projs
object as as_dict
method which offers a diret access to all the projected solutions. By default, it will enumerate the Boolean functions for each node. The method "count" instead returns the number of solutions per node. There is also a keys
parameter to specify a subset of nodes for the computation.
projs.as_dict(method="count")
{'Pax6': 1, 'Hes5': 1, 'Mash1': 1, 'Scl': 1, 'Olig2': 1, 'Stat3': 2, 'Zic1': 1, 'Brn2': 2, 'Tuj1': 13, 'Myt1L': 2, 'Sox8': 1, 'Aldh1L1': 1}
Note that the projected solutions gives an over-approximation of the full set of solutions: the full set of solutions is, in general, a strict subset of the cartesian product:
from functools import reduce
reduce(int.__mul__, _.values())
104
Access to the solutions of a specific node can be done as follows, where view
is an object similar to the one returned by bo.boolean_network()
(iterator over solutions).
with projs.view("Tuj1") as view:
functions = [f for f in view]
[str(f) for f in functions]
['Zic1', 'Myt1L|(Brn2&Zic1)', 'Zic1|(Brn2&Myt1L)', '(Brn2&Myt1L)|(Brn2&Zic1)', 'Myt1L|Zic1', 'Brn2&Zic1', '(Brn2&Myt1L)|(Brn2&Zic1)|(Myt1L&Zic1)', 'Brn2|Myt1L|Zic1', 'Brn2|(Myt1L&Zic1)', '(Brn2&Zic1)|(Myt1L&Zic1)', 'Brn2|Myt1L', 'Brn2', 'Brn2|Zic1']
Finally, projs
has a as_dataframe
method for pretty display of the projected solutions using pandas.
projs.as_dataframe()
Pax6 | Hes5 | Mash1 | Scl | Olig2 | Stat3 | Zic1 | Brn2 | Tuj1 | Myt1L | Sox8 | Aldh1L1 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Pax6 | !Mash1&Pax6 | !Hes5&Pax6 | Hes5&!Olig2 | Hes5&!Scl | Scl | Mash1 | Mash1&!Olig2 | (Brn2&Zic1)|(Myt1L&Zic1) | 0 | Olig2 | Stat3 |
1 | Hes5&Scl | Mash1 | Brn2|Myt1L | Olig2 | ||||||||
2 | Zic1 | |||||||||||
3 | Brn2 | |||||||||||
4 | Zic1|(Brn2&Myt1L) | |||||||||||
5 | Brn2|Myt1L|Zic1 | |||||||||||
6 | Brn2&Zic1 | |||||||||||
7 | Myt1L|Zic1 | |||||||||||
8 | (Brn2&Myt1L)|(Brn2&Zic1) | |||||||||||
9 | Myt1L|(Brn2&Zic1) | |||||||||||
10 | Brn2|Zic1 | |||||||||||
11 | (Brn2&Myt1L)|(Brn2&Zic1)|(Myt1L&Zic1) | |||||||||||
12 | Brn2|(Myt1L&Zic1) |
A self-contained bash-script for customizing the ASP code and performing the enumeration off-line can be exported as follows.
view = bo.boolean_networks()
view.standalone(output_filename="tutorial.asp")