Biodivine AEON#
This notebook illustrates how to use Python bindings of the tool AEON to symbolically analyze the attractors of non-trivial Boolean networks.
No special knowledge is necessary to understand the contents of this particular notebook. However, to successfully apply AEON to other problems, some basic understanding of AEON’s notion of partially specified Boolean network and the underlying symbolic representation using BDDs is recommended (these are explained here and here).
We start our exploration by loading a Boolean network from CellCollective:
from biodivine_aeon import *
import biodivine_aeon as ba
from pathlib import Path
import cellcollective
# Load the network from CellCellective.
sbml = cellcollective.load("https://research.cellcollective.org/#module/36604:1/signaling-pathway-for-butanol-production-in-clostridium-beijerinckii-nrrl-b598/10")
# Open the loaded SBML and parse it as a Boolean network.
model = BooleanNetwork.from_file(sbml.localfile)
#model = BooleanNetwork.from_file("cellcollective-36604-1.sbml")
print(f"Loaded model with {model.variable_count()} variables.")
Loaded model with 66 variables.
# Print some basic messages (unfortunately, in jupyter notebooks, these may be delayed until the cell finishes computing
# but you can see them in the terminal). ba.LOG_NOTHING is the default, but ba.LOG_VERBOSE is also available.
ba.LOG_LEVEL = ba.LOG_ESSENTIAL
Right now, the model has no explicit paramters, but there are several variables for which no update function is given in the SBML file. AEON will automatically consider these as implicit unknown parameters of the network:
print(f"Number of explicit parameters: {model.explicit_parameter_count()}.")
print(f"Number of implicit parameters: {model.implicit_parameter_count()}.")
Number of explicit parameters: 0.
Number of implicit parameters: 13.
Before we can explore the attractors of this network, first let us observe that the network as downloaded from CellCollective is actually not logically consistent with its regulatory graph. AEON will notify us about this fact when we try to construct the network’s asynchronous state-transition graph:
try:
stg = AsynchronousGraph(model)
except Exception as e:
print(e)
No update functions satisfy given constraints:
- spoIIE not activating in spoIIAB.
- sporulation has no effect in glucose___PTS.
- sigA has no effect in spo0A_p.
Fortunately, AEON can tell us exactly which properties are violated. So if we want to fix these issues, we know what variables to focus on. For now, we simply relax the conditions in the regulatory graph such that the existing update functions are valid:
# You can also automatically fix all issues using model.infer_valid_graph()
# Export the Boolean network as `.aeon` model string.
model_aeon = model.to_aeon()
# Relax the problematic regulation constraints:
# Mark regulation "spoIIE -> spoIIAB" as non-monotonous:
model.ensure_regulation({
'source': 'spoIIE',
'target': 'spoIIAB',
'essential': True,
'sign': None,
})
# Mark regulation "sporulation -| glucose___PTS" as non-essential:
model.ensure_regulation({
'source': 'sporulation',
'target': 'glucose___PTS',
'essential': False,
'sign': '-',
})
# Mark regulation "sigA -> spo0A_p" as non-essential:
model.ensure_regulation({
'source': 'sigA',
'target': 'spo0A_p',
'essential': False,
'sign': '+',
})
{'source': VariableId(51),
'target': VariableId(58),
'essential': True,
'sign': '+'}
Now we can actually compute the attractors quite easily. However, note that since AEON has to consider each possible variant of the network, the computation can still take several seconds or minutes (recall that we actually have 2^13 = 8192 possible networks).
stg = AsynchronousGraph(model)
attractors = Attractors.attractors(stg)
attractors
Start interleaved transition guided reduction with 604462909807314600000000[nodes:2] candidates.
[ColoredVertexSet(cardinality=867718363543552, colors=6704, vertices=867718363543552, symbolic_size=5924),
ColoredVertexSet(cardinality=15959465984, colors=1488, vertices=15959465984, symbolic_size=1393)]
> Discarded 302231454903657300000000 instances based on the BnVariable(63) extended component.
> State space reduced to 302231454903657300000000[nodes:5].
> Discarded 151115727451828650000000 instances based on the BnVariable(51) extended component.
> State space reduced to 151115727451828650000000[nodes:8].
> Discarded 75557863725914320000000 instances based on the BnVariable(55) extended component.
> State space reduced to 75557863725914320000000[nodes:10].
> Discarded 37778931862957160000000 instances based on the BnVariable(46) extended component.
> State space reduced to 37778931862957160000000[nodes:13].
> Discarded 18889465931478580000000 instances based on the BnVariable(45) extended component.
> State space reduced to 18889465931478580000000[nodes:16].
> Discarded 9444732965739290000000 instances based on the BnVariable(44) extended component.
> State space reduced to 9444732965739290000000[nodes:19].
> Discarded 4722366482869645000000 instances based on the BnVariable(43) extended component.
> State space reduced to 4722366482869645000000[nodes:22].
> Discarded 2361183241434822600000 instances based on the BnVariable(38) extended component.
> State space reduced to 2361183241434822600000[nodes:25].
> Discarded 1180591620717411300000 instances based on the BnVariable(37) extended component.
> State space reduced to 1180591620717411300000[nodes:28].
> Discarded 590295810358705700000 instances based on the BnVariable(35) extended component.
> State space reduced to 590295810358705700000[nodes:31].
> Discarded 295147905179352830000 instances based on the BnVariable(6) extended component.
> State space reduced to 295147905179352830000[nodes:34].
> Discarded 147573952589676410000 instances based on the BnVariable(5) extended component.
> State space reduced to 147573952589676410000[nodes:37].
> Discarded 73786976294838210000 instances based on the BnVariable(4) extended component.
> State space reduced to 73786976294838210000[nodes:40].
> Discarded 36893488147419103000 instances based on the BnVariable(3) extended component.
> State space reduced to 36893488147419103000[nodes:43].
> Discarded 18446744073709552000 instances based on the BnVariable(2) extended component.
> State space reduced to 18446744073709552000[nodes:45].
> Discarded 9223372036854776000 instances based on the BnVariable(36) extended component.
> State space reduced to 9223372036854776000[nodes:50].
> Discarded 4611686018427388000 instances based on the BnVariable(49) extended component.
> State space reduced to 4611686018427388000[nodes:61].
> Discarded 2448832297382707000 instances based on the BnVariable(59) extended component.
> State space reduced to 2162853721044680700[nodes:147].
> Discarded 937311672446484500 instances based on the BnVariable(60) extended component.
> State space reduced to 1225542048598196200[nodes:151].
> Discarded 468655836223242240 instances based on the BnVariable(62) extended component.
> State space reduced to 756886212374954000[nodes:155].
> Discarded 234327918111621120 instances based on the BnVariable(61) extended component.
> State space reduced to 522558294263332860[nodes:159].
> Discarded 117163959055810560 instances based on the BnVariable(53) extended component.
> State space reduced to 405394335207522300[nodes:167].
> Discarded 53409876830846980 instances based on the BnVariable(58) extended component.
> State space reduced to 351984458376675300[nodes:180].
> Discarded 18155135997837310 instances based on the BnVariable(64) extended component.
> State space reduced to 333829322378838000[nodes:204].
> Discarded 0 instances based on the BnVariable(25) extended component.
> Discarded 9288674231451648 instances based on the BnVariable(65) extended component.
> State space reduced to 324540648147386400[nodes:212].
> Discarded 18084767253659650 instances based on the BnVariable(52) extended component.
> State space reduced to 306455880893726700[nodes:217].
> Discarded 9042383626829824 instances based on the BnVariable(54) extended component.
> State space reduced to 297413497266896900[nodes:223].
> Discarded 4521191813414912 instances based on the BnVariable(56) extended component.
> State space reduced to 292892305453482000[nodes:229].
> Discarded 2260595906707456 instances based on the BnVariable(57) extended component.
> State space reduced to 290631709546774500[nodes:230].
> Discarded 356241767399424 instances based on the BnVariable(32) extended component.
> State space reduced to 290275467779375100[nodes:294].
> Discarded 178120883699712 instances based on the BnVariable(33) extended component.
> State space reduced to 290097346895675400[nodes:298].
> Discarded 72703339149656060 instances based on the BnVariable(39) extended component.
> State space reduced to 217394007746019330[nodes:841].
> Discarded 37279242387456 instances based on the BnVariable(29) extended component.
> State space reduced to 217356728503631870[nodes:1490].
> Discarded 1535450808320 instances based on the BnVariable(48) extended component.
> State space reduced to 217355193052823550[nodes:1496].
> Discarded 104155224670732290 instances based on the BnVariable(50) extended component.
> State space reduced to 113199968382091260[nodes:1494].
> Discarded 54338838528524290 instances based on the BnVariable(41) extended component.
> State space reduced to 58861129853566980[nodes:2031].
> Discarded 14433982238162944 instances based on the BnVariable(42) extended component.
> Finished ITGR process. 28 processes remaining.
> State space reduced to 44427147615404030[nodes:2095].
> Discarded 59408121856 instances using the BnVariable(13) transition basin.
> Finished ITGR process. 28 processes remaining.
> State space reduced to 44427088207282180[nodes:2218].
> Discarded 20658435540385790 instances based on the BnVariable(13) extended component.
> Finished ITGR process. 27 processes remaining.
> State space reduced to 23768652666896384[nodes:2630].
> Discarded 10329238934650880 instances based on the BnVariable(14) extended component.
> Finished ITGR process. 26 processes remaining.
> State space reduced to 13439413732245504[nodes:2670].
> Finished ITGR process. 26 processes remaining.
> Discarded 7746903680745472 instances based on the BnVariable(10) extended component.
> Finished ITGR process. 25 processes remaining.
> State space reduced to 5692510051500032[nodes:2431].
> Finished ITGR process. 25 processes remaining.
> Discarded 10468982784 instances based on the BnVariable(12) extended component.
> Finished ITGR process. 24 processes remaining.
> State space reduced to 5692499582517248[nodes:1678].
> Discarded 125829120 instances using the BnVariable(7) transition basin.
> State space reduced to 5692499456688128[nodes:1447].
> Discarded 1301624903434240 instances based on the BnVariable(7) extended component.
> State space reduced to 4390874553253888[nodes:2071].
> Finished ITGR process. 23 processes remaining.
> Finished ITGR process. 23 processes remaining.
> Discarded 309237645312 instances using the BnVariable(19) transition basin.
> Finished ITGR process. 23 processes remaining.
> State space reduced to 4390565315608576[nodes:2396].
> Discarded 3466340597760 instances based on the BnVariable(19) extended component.
> Finished ITGR process. 22 processes remaining.
> State space reduced to 4387098975010816[nodes:2525].
> Discarded 3270785368064 instances based on the BnVariable(34) extended component.
> Finished ITGR process. 21 processes remaining.
> State space reduced to 4383828189642752[nodes:2549].
> Discarded 1896420554833920 instances based on the BnVariable(27) extended component.
> Finished ITGR process. 20 processes remaining.
> State space reduced to 2487407634808832[nodes:4257].
> Finished ITGR process. 20 processes remaining.
> Discarded 645569491894272 instances based on the BnVariable(23) extended component.
> Finished ITGR process. 19 processes remaining.
> State space reduced to 1841838142914560[nodes:4393].
> Discarded 322784745947136 instances based on the BnVariable(28) extended component.
> Finished ITGR process. 18 processes remaining.
> State space reduced to 1519053396967424[nodes:4529].
> Discarded 319782912 instances based on the BnVariable(24) extended component.
> Finished ITGR process. 17 processes remaining.
> State space reduced to 1519053077184512[nodes:4655].
> Discarded 2916352 instances using the BnVariable(0) transition basin.
> Finished ITGR process. 17 processes remaining.
> State space reduced to 1519053074268160[nodes:4661].
> Discarded 46768661561344 instances based on the BnVariable(0) extended component.
> Finished ITGR process. 16 processes remaining.
> State space reduced to 1472284412706816[nodes:4527].
> Discarded 42135884922880 instances based on the BnVariable(30) extended component.
> Finished ITGR process. 15 processes remaining.
> State space reduced to 1430148527783936[nodes:3055].
> Discarded 7841709416448 instances using the BnVariable(11) transition basin.
> State space reduced to 1422306818367488[nodes:3491].
> Discarded 0 instances based on the BnVariable(11) extended component.
> Discarded 8126464 instances based on the BnVariable(40) extended component.
> State space reduced to 1422306810241024[nodes:3495].
> Discarded 6410457088 instances based on the BnVariable(21) extended component.
> State space reduced to 1422300399783936[nodes:3667].
> Discarded 3205228544 instances based on the BnVariable(17) extended component.
> State space reduced to 1422297194555392[nodes:3733].
> Discarded 1602614272 instances based on the BnVariable(18) extended component.
> State space reduced to 1422295591941120[nodes:3775].
> Discarded 801307136 instances based on the BnVariable(16) extended component.
> State space reduced to 1422294790633984[nodes:3817].
> Discarded 400653568 instances based on the BnVariable(1) extended component.
> State space reduced to 1422294389980416[nodes:3933].
> Discarded 865024 instances based on the BnVariable(47) extended component.
> State space reduced to 1422294389115392[nodes:3614].
> Discarded 370944 instances based on the BnVariable(22) extended component.
> State space reduced to 1422294388744448[nodes:3648].
> Discarded 240896 instances based on the BnVariable(31) extended component.
> State space reduced to 1422294388503552[nodes:3628].
> Discarded 129536 instances based on the BnVariable(20) extended component.
> State space reduced to 1422294388374016[nodes:3676].
> Discarded 64768 instances based on the BnVariable(15) extended component.
> State space reduced to 1422294388309248[nodes:3748].
> Discarded 44288 instances based on the BnVariable(9) extended component.
> State space reduced to 1422294388264960[nodes:3484].
> Discarded 4096 instances based on the BnVariable(26) extended component.
> State space reduced to 1422294388260864[nodes:3496].
> Discarded 2048 instances based on the BnVariable(8) extended component.
Interleaved transition guided reduction finished with 1422294388258816[nodes:3502] candidates.
Start Xie-Beerel attractor detection on 1422294388258816[nodes:3502] candidates.
> Found a bottom SCC: 867718363543552x6704[nodes:5924].
> Found a bottom SCC: 15959465984x1488[nodes:1393].
Attractor detection finished with 2 results.
As we can see, AEON found two separate attractor sets. Since the sets are very large, only a summary of each set is printed. Each set has millions of vertices and thousands of colors – AEON uses the term “color” to refer to concrete models arising for different parametrisations.
However, in this case, we can actually verify that these two sets are color-disjoint. This means that each possible network actually still only has one attractor. Consequently, we can just merge the two sets and treat them as a single attractor:
colors_in_both_sets = attractors[0].colors().intersect(attractors[1].colors())
print(f"Number of colors appearing in both sets: {colors_in_both_sets.cardinality()}")
attractor = attractors[0].union(attractors[1])
attractor
Number of colors appearing in both sets: 0
ColoredVertexSet(cardinality=867734323009536, colors=8192, vertices=867734323009536, symbolic_size=6154)
And just as we suspected, there is no intersection between the two sets. This can happen due to how the underlying algorithm works: some stable states are detected first, and more complex attractors are detected afterwards. However, since there are no intersections between these two sets, we can unify them and consider them as one attractor.
Now, to better understand what is happening in this attractor, we can classify the colors based on their long-term behavioural characteristics (this can also take a few second to compute):
classes = Classification.classify_long_term_behavior(stg, attractor)
classes
{Class(["disorder"]): ColorSet(cardinality=6144, symbolic_size=4),
Class(["stability"]): ColorSet(cardinality=2048, symbolic_size=4)}
In this case, we can see that 2048 colors correspond to a stable state, and 6144 colors correspond to a disordered attractor (i.e. an attractor that is neither a stable state nor a cycle). In this case, we have no oscillating (cyclic) attractors.
We can now use these classes to divide the attractor set based on its behaviour:
stable_attractor = attractor.intersect_colors(classes[Class("stability")])
disordered_attractor = attractor.intersect_colors(classes[Class("disorder")])
print(stable_attractor)
print(disordered_attractor)
ColoredVertexSet(cardinality=2048, symbolic_size=281)
ColoredVertexSet(cardinality=867734323007488, symbolic_size=5670)
To obtain some useful information from these sets, we can for example look at the possible values of specific network variables within each set. In this model, there is a variable called sporulation that is of particular interest, since it signifias that the cell is “hybernating”.
To study whether sporulation is on or off in the possible attractors, we can simply construct the set of all states where sporulation is on, and then compare it to the attractor sets:
sporulation_on = stg.mk_subspace({ "sporulation": True })
on_in_stable_attractor = stable_attractor.intersect(sporulation_on).vertices().cardinality()
off_in_stable_attractor = stable_attractor.minus(sporulation_on).vertices().cardinality()
print("Sporulation is ON in", on_in_stable_attractor, "stable states.")
print("Sporulation is OFF in", off_in_stable_attractor, "stable states.")
print("Sporulation is ON in", round((on_in_stable_attractor / (on_in_stable_attractor + off_in_stable_attractor)) * 100.0, 2), "% of stable attractor states.")
on_in_disorder_attractor = disordered_attractor.intersect(sporulation_on).vertices().cardinality()
off_in_disorder_attractor = disordered_attractor.minus(sporulation_on).vertices().cardinality()
print("Sporulation is ON in", on_in_disorder_attractor, "disorder states.")
print("Sporulation is OFF in", off_in_disorder_attractor, "disorder states.")
print("Sporulation is ON in", round((on_in_disorder_attractor / (on_in_disorder_attractor + off_in_disorder_attractor)) * 100.0, 2), "% of disordered attractor states.")
Sporulation is ON in 2048 stable states.
Sporulation is OFF in 0 stable states.
Sporulation is ON in 100.0 % of stable attractor states.
Sporulation is ON in 335223420223488 disorder states.
Sporulation is OFF in 532510902784000 disorder states.
Sporulation is ON in 38.63 % of disordered attractor states.
Interestingly, we have discovered that in every stable state, sporulation is active. We can thus use this fact to pick a witness network which guarantees that eventually, sporulation will always be active in a stable attractor.
This witness network is completely specified (i.e. has no parameters or uninterpreted functions). In other words, all of the explicit and implicit parameters are fixed to guarantee the desired behaviour:
# One color from the color set:
color_model = next(iter(classes[Class("stability")]))
witness = color_model.instantiate(model)
print(f"Number of explicit parameters: {witness.explicit_parameter_count()}.")
print(f"Number of implicit parameters: {witness.implicit_parameter_count()}.")
Number of explicit parameters: 0.
Number of implicit parameters: 0.
Finally, since computing the attractors can take a long time, it may be necessary to save them into a file for further processing. For this, we can export the underlying BDD representation:
Path("stable_attractor.bdd").write_text(stable_attractor.to_bdd().data_string())
Path("disordered_attractor.bdd").write_text(disordered_attractor.to_bdd().data_string())
62415
Or we can use pickle:
import pickle
with open("stable_attractors.bdd.pickle", "wb") as f:
pickle.dump(stable_attractor.to_bdd(), f)
with open("disordered_attractor.bdd.pickle", "wb") as f:
pickle.dump(disordered_attractor.to_bdd(), f)
With the BDDs saved, we can always reload them from disk. However, remeber that the BDD does not contain any information about the Boolean network from which it was constructed! So you must make sure that you only import BDDs that were created using the same Boolean model. Ideally, always save an .aeon or .sbml file along with any raw BDD data.
network_context = stg.symbolic_context()
bdd_context = network_context.bdd_variable_set()
stable_reloaded = ColoredVertexSet(network_context, Bdd(bdd_context, Path("stable_attractor.bdd").read_text()))
disordered_reloaded = ColoredVertexSet(network_context, Bdd(bdd_context, Path("disordered_attractor.bdd").read_text()))
print(stable_reloaded)
print(disordered_reloaded)
assert stable_reloaded == stable_attractor
assert disordered_reloaded == disordered_attractor
ColoredVertexSet(cardinality=2048, symbolic_size=281)
ColoredVertexSet(cardinality=867734323007488, symbolic_size=5670)
As we can see, the sets are the same as the ones that we saved into the files.
Finally, let us quickly note that using AEON, you can also compute reachability from/to a specific symbolic set. Furthermore, you can specify a subset within which the reachability should be computed. Using this mechanism, you can often build more complex analysis techniques. For example, you can analyse the strong and weak basin of a particular attractor, or check for reachable attractors from a specific initial state.
# Pick a set containing a single arbitrary vertex from the whole state space:
vertex = stg.mk_unit_colored_vertices().pick_vertex()
# The actual vertex is just a list of Boolean values.
print("Vertex", next(iter(vertex.vertices())).values())
# Compute the set of states backward-reachable from the vertex.
basin = Reachability.reach_bwd(stg, vertex)
print("Basin", basin)
# Compute the strongly-connected-component in which the vertex resides.
scc = Reachability.reach_fwd(stg, vertex).intersect(basin)
print("SCC", scc)
Vertex [False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False]
Basin ColoredVertexSet(cardinality=11558771119839541985280, symbolic_size=101)
Start basic symbolic reachability with 8192[nodes:68] candidates.
SCC ColoredVertexSet(cardinality=8192, symbolic_size=68)
Basic reachability done: 11558771119839542000000[nodes:101] candidates. Max intermediate size: 441.
Start basic symbolic reachability with 8192[nodes:68] candidates.
Basic reachability done: 44364707639821350000[nodes:28057] candidates. Max intermediate size: 43163.
For example, for an initial state where all variables are inactive, we can see that its (weak) basin is quite large, but the vertex is in fact a trivial SCC (the whole strongly connected component consists only of one vertex, regardless of color).
So far, this is the end of our demo. In case of any issues, feel free to contact us on github!
To learn more about what features and functions are available in AEON.py, you can follow the tuturial on partially specified Boolean networks and symbolic computation using BDDs. Further tutorials, examples and case studies are available in the AEON.py repository. API documentation is also available here.